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We study, analytically and with lattice simulations, the decay of coherent field oscillations and the 
subsequent thermalization of the resulting stochastic classical wave-field. The problem of reheating 
of the Universe after infiation constitutes our prime motivation and application of the results. We 
identify three different stages of these processes. During the initial stage of "parametric resonance" , 
only a small fraction of the initial infiaton energy is transferred to fluctuations in the physically 
relevant case of sufficiently large couplings. A major fraction is transfered in the prompt regime 
of driven turbulence. The subsequent long stage of thermalization classifies as free turbulence. 
During the turbulent stages, the evolution of particle distribution functions is self-similar. We 
show that wave kinetic theory successfully describes the late stages of our lattice calculation. Our 
analytical results are general and give estimates of reheating time and temperature in terms of 
coupling constants and initial infiaton amplitude. 
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I. INTRODUCTION 

Field theoretical systems which are a long way from 
thermal equilibrium have been studied intensely in re- 
cent years. The particular problem of how and when 
such systems approach equilibrium stretches beyond ob- 
vious fundamental interest and finds many practical ap- 
plications. In high-energy physics understanding of these 
processes is crucial for applications to heavy ion collisions 
and to cosmology of the early universe. The first topic 
gains further importance in light of the current and future 
experimental search for a quark-gluon-plasma at RHIC 
and at the forthcoming LHC. The second application, our 
main interest in this paper, is related to the problem of 
reheating of the universe after cosmological inflation. 

Inflation provides a solution to the flatness and the 
horizon problems of standard cosmology P, S 13 and 
explains the generation of initial density perturbations 
- the seeds of galaxies and large-scale structure in our 
universe. During inflation the universe is in a vacuum- 
like state. At the end of inflation all energy density is 
stored in a Bose condensate, the coherently oscillating 
"infiaton" field. This state is highly unstable: paramet- 
ric, tachyonic or strong non-adiabatic particle creation 
triggers a fast and explosive decay of the infiaton. This 
process, dubbed preheatin g H , | a , is currently well un- 
derstood [iSSIlEOillllll- A generic feature is a 
strong and fast amplification of fluctuation fields at low 
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momenta, which may lead to various interesting phys- 
ical effects during and after preheating. These include 
non-thermal phase transitions [T^ ITsL ll^ ll^ w ith pos- 
sible formation of topological defects Il8lll9ll20ll2ll l2^. 
creation of super- heavy particles [231. generation of 
high-frequency gravitational waves p4| . etc. 

The explosive stage of infiaton decay ends when the 
rate of interactions of created fluctuations among them- 
selves and with the infiaton becomes comparable to the 
infiaton decay rate |E S IE 0- The understanding of 
the subsequent stages of relaxation towards equilibrium, 
of the thermalization processes and the calculation of 
the final equilibrium temperature is important for var- 
ious applications as it links the inflationary phase with 
that of standard cosmology. Among those one can list 
baryogenesis (1^ .2^ 23t I^El and the problem of over- 
abundant gray itino production in supergravity models 
|2llMl3i r l3 l l3ll33|. It determines the abundances of 
other relics, like super- heavy dark matter (35- .36.. .37. ..38.J . 
or axino dark matter [39l |. 

Knowledge of the reheating temperature is also im- 
portant for fixing constraints on the inflationary model 
from Cosmic Microwave Background (CMBR) anisotropy 
|40ll4lll4^ . In some models cosmologically important 
curvature perturbations may be even generated during 
the process of thermalization 0, 0, [4^ @ • Last 
but not least: the reheating temperature should be larger 
than about one GeV to ensure that the standard Big 
Bang Nucleosynthesis 0;E3 not hampered. 

There have been many efforts and successes in the un- 
derstanding of the non-equilibrium dynamics and relax- 
ation of field theories, see e. g. Refs . [sol IsiL Is^ IsM 
iS El Isl El lell^ However, 
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the leading asymptotic dynamics towards equilibrium re- 
mained rather less understood and developed. 

The main problem for the theoretical understanding 
of reheating is that initially the occupation numbers are 
very large, of order of the inverse coupling constant. In 
addition, in many inflationary models the zero mode does 
not decay completely during preheating. Therefore, a 
simple perturbative approach is not justified. On the 
other hand, in this regimeja description in terms of clas- 
sical field theory is valid and the whole process (in- 
cluding preheating), can be studied by classical lattice 
simulations. 

Recently, we employed this method to show [g^ 
that the classical reheating of a massless $^-theory in 
3-1-1 dimensions is characterized by a turbulent and self- 
similar evolution of distribution functions towards equi- 
librium. The shape of the spectra, as well as the self- 
similar dynamics, could be understood within the frame- 
work of wave kinetic theory. This made it possible to 
estimate reheating time and temperature, which turned 
out to coincide parametrically with the results of the sim- 
ple perturbative approach. 

Turbulence appears in a large variety of non- 
equilibrium-phenomena in nature (see Refs. |68l l69l l70l | 
for a general introduction). It was first discussed for flu- 
ids, in the regime of large Reynolds numbers (velocities), 
where viscosity is subdominant. Kolmogorov identified 
turbulence in this regime |^ [t^ as a statistically scale 
invariant flow of spectral energy mediated by vortex in- 
teractions. The same dynamical structure may appear 
in systems of coupled wave s, e . g. o n fluid surfaces or for 
coupled flelds in a plasma l69l 173^. In this case the 
cascade is mediated by wave interactions and the phe- 
nomenon has been called wave turbulence. 

If there exists an active (stationary) source of energy 
in momentum space, the turbulence is called driven (sta- 
tionary) . When the source is switched off after the stage 
of activity, the freely propagating energy cascade is often 
referred to as free turbulence. If the kinetic description is 
applicable, the energy cascade is called weak turbulence. 
Otherwise one is facing a strong turbulence. 

One may expect that the concept of turbulence should 
be relevant for the problem of reheating |^[53| already on 
general grounds. Indeed, the source of energy, localized 
in the "infra-red" is present initially. It is represented 
by the inflaton held in the problem at hands. To com- 
plete the argument, we note that as the flnal outcome 
of the evolution one should expect cascading of energy 
towards a signiflcantly separated region of "ultra-violet" , 
high momentum modes. 

The goal of the current paper is twofold. First, we 
want to apply the wave kinetic theory of turbulence to 
the problem of Universe reheating after inflation. We 
derive general formulas for the spectra of turbulent dis- 
tributions and for the self-similar evolution towards equi- 
librium. This enables us to give asymptotic estimates of 
reheating time and temperature in Minkowski space as 
well as in Friedmann universe. 



Second, we want to test and confront these ideas to 
numerical lattice calculations. For our numerical integra- 
tions we have chosen the simplest "chaotic" inflationary 
model [t^I • While the initial "preheating" stage in other 
inflationary models, e.g. in hybrid inflation |75| may ex- 
hibit important differences IzS with this model, 
we expect the subsequent turbulent stages to be more 
universal. 

We start lattice integration from "vacuum" initial con- 
ditions for fluctuations in a background of classical oscil- 
lating inflaton field. We observe the initial parametric 
resonance stage when the energy in fluctuations is grow- 
ing exponentially with time. This stage terminates when 
re-scattering of waves out of the resonance band becomes 
important '3 • physically relevant case of suf- 

ficiently large couplings this happens rather early, when 
only a small fraction of initial inflaton energy is trans- 
ferred to fluctuations 0, 0, . At this point a state of 
stationary turbulence should be established that is driven 
by the zero- mode. On general grounds, it can be deduced 
that during this stage the energy in fluctuations should 
grow linearly with time. This behavior is conflrmed by 
the results of our numerical simulations. The stage of 
stationary turbulence should terminate when the energy 
left out in the zero-mode becomes smaller than the en- 
ergy stored in created " particles" . From this moment of 
time, the transport of energy from the source is negligible 
and we observe free turbulence with self-similar evolution 
of particle distributions towards thermal equilibrium. 

The first stage of driven turbulence is prompt and gives 
the main mechanism by which energy is drawn out of the 
zero-mode, e.g. out of the inflaton fleld. The identifi- 
cation of this constitutes one of the new results of the 
present paper, as opposed to the common opinion that 
the main mechanism is a " parametric resonance" . The 
second stage of free turbulence is very long and can be 
analytically described as self-similar evolution. This is 
another new result and diffuses some existing claims and 
hopes that "parametric resonance" may bring a system 
to thermal equilibrium on a very short time scale. 

Overall, the kinetic description and the results of lat- 
tice simulations are in rather good agreement with each 
other. This indicates that the regime of weak wave tur- 
bulence may be already achieved on the lattice. 

The paper is organized as follows. In Section^we re- 
view the results of our numerical simulation of reheating 
in the simplest model to get familiar with concepts, 
problems and the typical dynamical behavior of the sys- 
tems of interest. In Sections IIIII IIVI we apply the theory 
of wave turbulence to the problem of reheating in gen- 
eral. In Section we present our numerical simulations. 
In Section IVTI we compare lattice results with the kinetic 
approach and discuss the applicability of the latter. In 
Section IVIll we discuss some physical applications of our 
results, in particular the thermalization in the self-similar 
regime. In Appendix 1X1 we give the details of our numer- 
ical procedure. In Appendix IbI we review the derivation 
of the kinetic equation for a system of weakly interacting 
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classical waves. 



II. THE SIMPLEST MODEL OF REHEATING. 
NUMERICAL RESULTS. 

We start with a presentation of our numerical results 
for the inflaton decay and the subsequent equilibration of 
the decay products in a simple A$'* model. The results 
were already briefly reported in Ref . 0| . The numerical 
procedure itself is described in Appendix IXI At the end 
of the Section we will discuss some expected differences 
with more complicated models. This order of presenta- 
tion allows us to introduce the typical behavior in the 
systems under consideration and to formulate concepts 
and problems. This will be useful in the discussion of 
the general theory of turbulent thermalization, which we 
carry out in the following Section. Further numerical 
results, obtained for the simplest model, and nu- 
merical results obtained for more complicated multi-field 
systems, will be presented in Sections Ivl and IVTl 

A. Results for the $^-Model 

In this simple model, the field $ is the only dynamical 
variable. Its initial homogeneous mode drives inflation, 
while development and growth of fluctuations on sub- 
horizon scales at the end of inflation can be viewed as a 
simple model of reheating. Inflation ends when the mo- 
tion of the homogeneous mode of the field changes from 
the regime of "slow-roll" to the regime of oscillations. It 
is convenient to work in conformal coordinates where the 
metric takes the form 

ds^ =a{'qf{drf -dx^). (1) 

We choose the case of a massless field where the equation 
of motion for the rescaled field Lp = $a after inflation is 
the same as in flat space-time 

Dip + Xip^ ^0. (2) 

Therefore, all results obtained in this model are equally 
applicable to the reheating of the Universe after inflation 
and to modeling of other processes of thermalization in 
relativistic systems, say, after heavy ion collisions. 

The homogeneous component of the field, which corre- 
sponds to the zero momentum in the Fourier decomposi- 
tion, ipo{ri) = {(f), is usually referred to as the "zero- 
mode." It is convenient to make a rescaling of the 
field, = ip/ipo{r]Q), and of the space-time coordinates, 
x'^ ^/\(po{r]o)xf^ . Here, t]o corresponds to the initial 
moment of time (end of inflation). In what follows di- 
mensionless time is still denoted as rj. With this rescal- 
ing, the initial condition for the zero-mode oscillations is 
00 ('yo) = Ij 3,nd the equation of motion takes the simple 
parameter free form 

□0 + ,/,3^O. (3) 
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FIG. 1: Squared amplitude of the zero-mode oscillations, (/!)q, 
and variance of the field fluctuations as functions of time 77. 



All model dependence on the coupling constant A and 
on the initial amplitude of the field oscillations is en- 
coded now in the initial conditions for the small (vac- 
uum) fluctuations of the field with non-zero momenta 
(see and appendix ^ . The physical normalization 
of the inflationary model corresponds to a dimensionful 
initial amplitude of (/3o('7o) ~ 0.3Mpi and a coupling con- 
stant A ~ 10"^'^ The re-parameterization property 
of the system allows to choose a larger value of A for 
numerical simulations. We have used A = 10^*. 

Various quantities can be measured in a lattice cal- 
culation and monitored as functions of time. Here we 
will discuss the zero mode, 0o = the variance, 

var((j)) = {(jy^) — ipQ and "particle occupation" numbers. 
For definitions see Appendix fXl 

We begin the discussion of our numerical results with 
the evolution of the zero-mode and the variance of the 
field, which are shown in Fig. ^ The zero mode is a 
rapidly oscillating function on the time scale of our lat- 
tice calculation. In Fig. we show the amplitude of 
oscillations, 0g, as a function of time. 

The initial fast transfer of the zero-mode energy into 
fluctuations during preheating (up to 77 ~ 300) is fol- 
lowed by a long and slow relaxation phase. In this late 
time regime the amplitude of the zero mode oscillations 
decreases according to ^ rj~^ with z ~ 1/3, the variance 
of the field (averaged over high-frequency oscillations) 
drops according to ~ rj~^ with v « 2/5. In addition, we 
find that in this regime the zero-mode is in a non-trivial 
dynamical equilibrium with the bath of highly occupied 
modes: when the zero-mode is artificially removed, it is 
recreated on a short time-scale (Bose condensation). 

A detailed analytical discussion of the initial linear 
stage of the parametric resonance in this model can be 
found e.g. in Refs. [tI, |7i,|73. During this stage the 
occupation numbers grow exponentially with time in a 
narrow band of resonance momenta. Figure [21 shows the 
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FIG. 2: Occupation numbers as function of k/cpQ at 
ri = 100, 400, 2500, 5000, 10000. 



occupation numbers at different moments of time. The 
displayed spectrum at time r] = 100 corresponds to the 
stage of parametric resonance. The resonance peak is lo- 
cated at the theoretically predicted value of fcrcs. ~ 1.27 
ff^ . At later time, the growth of the resonance peak is 
stopped by re-scattering of particles out of the resonance 
band, which leads to a broadening of the occupied region 
and to the appearance of multiple peaks of compara- 
ble width, see spectrum at ?7 = 400. This structure fits 
estimates for the development of turbulence in the pres- 
ence of a narrow width source located at a finite k, see 
Ref. |6^. At even later times the spectra have become 
smooth because of re-scattering, and only the first peak 
is still visible as a small bump. With time, its position 
moves towards smaller momenta, reflecting the change in 
the effective frequency of inflaton oscillations. However, 
if the particle momenta are rescaled by the current am- 
plitude (j)o of the zero mode oscillations, as in Fig. |21 the 
position of the resonance is approximately unchanged. 
Particles with small momenta are distributed according 
to a power law, which at larger momenta is bounded by 
a cut-off. The position of this cut-off moves with time 
to the " ultra-violet" . This reflects a general tendency of 
the system to thermal equilibrium. Indeed, in a state of 
thermal equilibrium the energy of the system should be 
concentrated at much higher wave-numbers compared to 
the resonance momenta. On the other hand, energy is 
inputted into the system of particles in the region of k 
near the resonance peak. Therefore, we have a continu- 
ous flux of energy across momentum space, from low to 
high momenta. 

This stage of evolution (77 > 1500) has the following 
characteristic features : 

1. The system overall is statistically close to a Gaus- 
sian distribution of field amplitudes and conjugated 
momenta l57ll6ql. 



FIG. 3: On the right hand side we plot the wave en- 
ergy per decade found in lattice integration at 77 = 
3600, 5100, 7000, 10000. On the left hand side are the same 
graphs transformed according to the relation inverse to 
Eq. @. 



2. The spectra in the dynamically important region 
can be described by a power law, fc^'* with s « 3/2. 
We see that the system is not in a thermal equilib- 
rium which would correspond to s = 1. Rather, the 
exponent of particle distributions in the power law 
region corresponds to Kolmogorov turbulence 66]. 

3. The power law is followed by a cut-off at higher k. 
Energy accumulated in particles is concentrated in 
the region were the cut-off starts. Its position is 
monotonously growing toward the ultra-violet, re- 
flecting the evolution towards thermal equilibrium. 

4. This motion can be described as a self-similar evo- 
lution ii 



n{k,T) = T-^noikT-P) , 



(4) 



where r = rj/rjc and rjc is some (arbitrary) late-time 
moment. The best numerical fit corresponds to g « 
3.5p and p w 1/5, and is presented in Fig. O The 
value of the exponent p is of prime interest since it 
determines the rate with which system approaches 
equilibrium. 

The first and the second point in this list facilitate the 
use of wave kinetic theory, see e.g. |6^|83. However, a 
straightforward application is difficult and may be even 
inappropriate, at least at the early re-scattering stages, 
because of the following: 

1. The zero mode does not decay completely. It may 
induce "anomalous" terms in the collision integral, 
which are absent in the usual kinetic description. 

2. The occupation numbers are large initially, of or- 
der of the inverse coupling constant, Uk l/X, see 



5 



Fig. 121 Therefore, in addition to fowest order colli- 
sions (e.g. scattering of two particles into two par- 
ticles with different momenta), multi-particle colli- 
sions may be dynamically important as well. 

Therefore, precise lattice calculations are needed. On 
the other hand, they have a limited dynamical range in 
momenta and in time, and one has to switch to a kinetic 
approach at some later stage. To determine when (and 
if) this is possible, the results obtained with the use of 
a simple kinetic approach should be confronted with the 
lattice results. 

In the following Sections we will develop and apply 
the theory of weak wave turbulence to the models of the 
type integrated on the lattice. In particular, we will cal- 
culate all universal scaling exponents and show that they 
are in agreement with lattice results. At "early" times 
the dynamics of the model described above is driven by 
m-particle interactions with to = 3. Wave turbulence 
theory gives for scaling exponents in d = 3 spatial di- 
mensions: 

p= l/(2m- 1) , 

s = d — m/{m — 1) , 

V = 2/(2m- 1) , 

z = 1/ {d{m — 1) — to) . 

B. Expected differences in more complicated 
models 

The flux of energy over momentum space, which is nec- 
essarily present in problems like reheating and thermal- 
ization after inflation, signifles that we should observe a 
turbulent state during the thermalization stage and that 
the theory of turbulence applies. In a simple model 
the stage of preheating (i.e. parametric resonance) ends 
when roughly half of the inflaton energy is transferred 
to particles. Indeed, Fig. ^ shows that the amplitude of 
the zero mode, which is a source of energy for the tur- 
bulence problem, starts to decrease already at the end of 
the parametric resonance stage. In such system we ex- 
pect the free turbulence regime to follow the preheating 
stage. 

In more complicated systems, which involve other 
fields coupled to the inflaton, say, some field Xi paramet- 
ric resonance may end when the fraction of energy trans- 
ferred to the X excitations is still negligible compared 
to the energy stored in the infiaton zero mode. Indeed, 
parametric resonance ends when the rate of re-scattering 
of particles out of the resonance band became compa- 
rable to the resonant production rate and the maximal 
value of the variance of x excitations achieved at the end 
of the resonance stage is ~ l/ff^, where is either the 
coupling of X to the infiaton, or self-coupling of x (viz., 
the largest of these two). We expect that in this case 
turbulent transport will develop when the amplitude of 
the inflaton zero mode is still unchanging. This means. 



that the transfer of zero mode energy into ^-filed should 
occur in the regime of stationary turbulence. Only when 
the amount of energy in the zero- mode becomes subdom- 
inant we should expect a transition to the regime of free 
turbulence. This is an important difference to the sim- 
ple (/)^-model. In particular, the distribution functions 
move much faster into the ultra-violet in this regime, 
p = (to — 1)/(2to— 1). We will see that the regime of sta- 
tionary turbulence is indeed present in two field models, 
see Section [3 

III. THERMALIZATION IN THE WAVE 
KINETIC REGIME. GENERAL THEORY. 

A. Turbulent reheating: a motivation 

Kolmogorov's turbulence is characterized by a station- 
ary transport of some conserved quantity between differ- 
ent scales in momentum (Fourier) space |7ll l72j . In the 
following, we will restrict ourselves to systems with spa- 
tially isotropic and homogeneous correlation functions, 
which applies to the cosmological conditions after infia- 
tion. Turbulence usually appears when a source of energy 
or particles is present and is localized in some momen- 
tum region /cjn. In addition to the source exists a "sink" 
localized at /cout- When both, source and sink are station- 
ary, it is natural to expect the eventual development of a 
stationary state with scale independent transport of the 
conserved quantity through momentum space. Indeed, 
energy or particle number cannot accumulate between 
/cin and kont and should fiow from one scale to the other. 

This is a system-independent formulation of Kol- 
mogorov's concept of turbulence, which he formulated 
in the context of hydrodynamical systems ItiI l73l . Za- 
kharov applied it to systems of coupled waves |73 in the 
regime of kinetic wave interactions. His approach is based 
on his derivation of the wave kinetic equations (see e.g. 
[iElzEliOl) and is well suited to studies of turbulence in 
classical field theories. We will adopt it here. 

The physical scenario of reheating after infiation shares 
basic ingredients with that of turbulence: there exists a 
localized source of energy- the coherently oscillating in- 
flaton zero-mode - pumping energy into the system of 
particles at Fourier wave- numbers ki^ ~ fcrcs- The mech- 
anism behind this pumping can be parametric resonance, 
tachyonic amplification, etc. Like in the turbulent sce- 
nario there do not exist other intermediate scales (wave- 
numbers), where energy is infused, accumulated, or dissi- 
pated. Thus, it seems likely that the eventual dynamics 
of reheating - after the explosive regime of preheating has 
ended - is close to that of Kolmogorov's turbulence. 

However, in the description of reheating appear some 
differences to stationary turbulence, since: 

1. A sink does not exists. 

2. The source (i.e. the amplitude of inflaton zero- 
mode oscillations and therefore interaction rates) 
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can be essentially time-dependent on relevant time 
scales. 

3. Neither source nor sink exist when the inflaton has 
completely decayed. 

In the first case, we expect that the stationary turbu- 
lent flux of energy still will be established in some "iner- 
tial" range fcin < k < fcout- Particle distributions in this 
range of momenta should not be significantly different 
compared to the case with a stationary sink. Indeed, in 
the typical turbulent problem the energy dissipates (e.g. 
into heat) after entering the region k 3> fcout- For prob- 
lems relevant to thermalization after inflation, instead of 
dissipation the transported energy is used to populate 
high momentum modes at fc 3> fcout- If the transport is 
reasonably "local" in momentum space, the flux of energy 
through the inertial range should not be influenced much 
by processes which involve fc > fcout- Energy may dissi- 
pate at fcout or continue the flow to even higher momenta, 
but regardless of this, we should expect the same distri- 
bution of particles in the inertial range. However, in the 
latter case we can expect that the value of fcout increases, 
and since the flux of energy is constant throughout in- 
ertial range, the total energy of a system without a sink 
has to grow linearly with time. 



E{t) cx t . 



(5) 



This is a simple consequence of the stationarity of tur- 
bulence in the inertial range, and can be used as its sig- 
nature. 

A time dependent source (second point above) changes 
the picture somewhat, since stationary states are not 
likely to develop even in a finite range of fc. However, 
a weak time-dependence should still allow for a close-to- 
stationary and close-to-turbulent evolution. Moreover, 
even if the source eventually does not exists, particle dis- 
tributions in the inertial range as functions of momenta 
can still be close to turbulent power laws. Indeed, sta- 
tionary turbulent distributions can be found as zeros of 
the collision integral - In the non-stationary case the 
collision integral is non-zero, but should approach a min- 
imal value in the inertial range which may result in the 
same shape of particle distributions there. 



B. Wave turbulence by scaling analysis 

The dynamics of coupled waves close to a stationary 
state can be described by a wave kinetic equation (see 
e.g. refs. [HItIIP): 



fik = h [n] 



(6) 



Here the function nk, usually called occupation number 
or wave action, describes the average volume of phase 
space occupied by the oscillations of a single mode with 
a wave-number fc. Its evolution is a result of resonant 
wave interactions, the effect of which is described by the 



collision integral Ik [n] . The collision integral is a function 
of the "external" momentum fc and a functional of the 
distribution function n, which is reflected in the notations 
we use. When we do not need to stress the functional 
dependence, we will also write Ik as /(fc). The collision 
integral for the case of interest, Eq. jSJ, is explicitly 
derived in Appendix IbI 

Before we proceed, let us remind the general structure 
of the collision integrals using as illustration the scat- 
tering of two particles into two particles, which will be 
referred to as 4-particle process. This will also allow us 
to introduce the necessary notations. In all cases we will 
write the collision integral as 



Ik[n] - / dn{k,q,)F{k,q,). 



(7) 



This form separates the contributions which are due to 
the (fixed) particle model, dfl{k,qi), from those which 
are due to the (evolving) particle distribution functions, 
F(k,qi). Here fc is the external momentum and qt refer to 
momenta over which the integration is carried out. If m 
particles participate in the collision, i takes values from 
1 to TO — 1. E.g. when 2 particles scatter into 2 particles, 
m — A and there are 3 internal momenta over which we 
integrate, qi, q2 and q^. Namely 



(27r)^|Af|% ,^,^ ' d% 

dil[k,qi)^ — d [kf,,qif,) [[■ 



2ujk 



-I J- 2w,(27r)3 

Z— 1 ^ ' 



(8) 



dn contains the usual energy-momentum conservation 5- 
functions, which we have denoted as d*{k^, lifi), the "ma- 
trix element" squared, |Mp, of the corresponding process 
(which is a function of fc and qi) and the integration mea- 
sure over momentum space. Here, fcg = Wfc = CL'(fc) and 
uji = uj{qi) refers to the particle energy. 

When quantum effects are accounted for, the function 
F in our example is given by 

F{k,qi) = (1 + Uk) (1 + nq^)nq^nq^ 

- nkUqJl + Uq^) {1 + riq^) . (9) 

In the limit n ^ 1 terms 0{n^) can be neglected and F 
is a sum of terms 0{n^) 

F{k,qi) = {uk +nq^)nq^nq^ - nknq^{nq^ +^^93) - (10) 

The limit n ^ 1 corresponds to interaction of classical 
waves and expression Eq. H10|) is also explicitly derived in 
Appendix FBI This illustrates a general rule: in the classi- 
cal limit and for interaction of m waves the function F is 
a sum of terms 0(n"^~^) with appropriate permutations 
of signs and indices. In other words, in this limit _F is a 
homogeneous function with respect to multiplication of 
each occupation number by Q 



F{Cn) = r-'F(Cn) 



(11) 



This property is extremely important in our subsequent 
analysis. When quantum effects become important (i.e. 
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when one should properly write [1 + n] in F), the clas- 
sical turbulence and/or self-similar evolution stops. At 
that moment particle distributions relax to usual Bose- 
Einstein functions. We will not be concerned here with 
the (presumably relatively short) relaxation period from 
the classical to the quantum regime, but will study in 
detail the turbulent evolution in the regime of classical 
waves. 

This gives us sufRcient notational details to proceed 
with the discussion of turbulence. We restrict it to sys- 
tems which are isotropic and homogeneous in configu- 
ration space, when occupation numbers (as well as all 
other parameters which enter the collision integral) de- 
pend on the modulus of momenta only. We consider the 
classical limit in the function F with general m-particle 
interaction, in case of which Eq. Hll|) holds. To keep 
the discussion general, in the rest of this section we will 
consider the case of (d-l-1) dimensional space time. 

Often a collision integral conserves one or several quan- 
tities. We restrict ourselves to energy density 



P 



{2tt 



(12) 



which is conserved when the expansion of the Universe 
can be neglected or "rotated" away, and particle density 



(2^ 



(13) 



which corresponds to conserved charges, e.g. baryon 
number. 

Conservation of n or p can be expressed as a continuity 
equation in Fourier space, e.g. 



dt{ujk nk) + Vfe • jk 0. 



(14) 



Here and in what follows we will write the explicit re- 
lation for energy conservation, the case of conserved 
charges can be easily obtained by a formal substitution 
Wfc = 1. In the isotropic case only the radial component 
of the flux density, jk, is non- vanishing and we get for 
the energy flux, S''{p), trough the sphere of radius p 



r(i + f) 



(15) 

dkk'^ ^ujkhin] , 



In (|15|l the factor in front of the integral is the area of the 
d-dimensional unit sphere. In case of stationary turbu- 
lence this flux should be scale-independent, i.e. integral 
Eq. p5|l should not depend upon its integration limit p. 
This is possible if the collision integral equals zero. One 
can exphcitly look for solutions Ik[n] — 0, see e.g. [6^ . 
Such solutions correspond to stationary turbulence and 
exist with non-trivial boundary conditions (source and 
sink), in addition to the Rayleigh- Jeans-law of classical 



equilibrium. Here we adopt an alternative and somewhat 
simpler approach of Ref. 82] to determine the turbulent 
solutions. 

Following 82] we consider states for which the collision 
integral has certain scaling properties under ^-rescaling 
of the external momentum k 



(16) 



To simplify notations we assume that all momenta were 
made dimensionless by rescaling with some typical mo- 
mentum scale, without explicitly writing this. The spe- 
cial choice ^ = fc~^ allows us to find the fc-dependence 
of the collision integral, Ik[n] = A:~'^/i[n]. Let us addi- 
tionally assume that the dispersion law is a homogeneous 
function as well, 



(17) 



Relations 1)1 6(1 and H17|l should hold in some region of mo- 
menta where we expect turbulent behavior. Integrating 
Eq. we find 



S{p) (X --'^+°-'^ 



P 



d + a — 



(18) 



Here we indicated explicitly that the collision integral 
in the turbulent state with scaling behavior Eq. (|16|) 
depends on the exponent v. We find that the flux is 
scale invariant, if 



v ^ d + a . 



(19) 



This condition defines the turbulent exponents which we 
will specify in detail below. Note that this implies the 
existence of the limit 



lim 

u — >d-\-a 



const 7^ 



(20) 



as a sufhcient condition for the existence of a stationary 
turbulent solution: if the collision integral has a zero of 
first degree at z/ = d -I- a , the turbulent flux is scale- 
invariant and finite. 

In what follows, we consider particle models for which 
dfj is a homogeneous function of all momenta 

dn{^k,iq,)=edn{k,q,) . (21) 

Rescaling of the external momentum fc by ^ gives 



kk^e j d^{k,q,)F{ik,iq,), 



(22) 



since integration over every qi is from to cx3. We will 
exploit this relation in two ways: 

1. Often the evolution of distribution functions in- 
volves rescaling of their momenta, see Sec. 1111 CI 
If this is the case, the collision integral as a func- 
tion of time can be found with the help of 



dn{k,q,)F(ikXq,)^C''hk 



(23) 
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2. Let us assume that the particle distribution func- 
tions are power laws in the momenta, 

n{q) cx q-' . (24) 

This leads to the following scaling of F 

FiCk,iq,)^r'("^''^F{k,q,), (25) 

Combining this with Eq. ^ we find I^k = C^"" h- 
A comparison with Eqs. and p9|l leads us to the ex- 
ponent s which defines the scaling of particle distribution 
functions in a turbulent state with constant energy trans- 
port (we will call this energy cascade for brevity) 



(26) 



m — 1 



Turbulence with constant transport of particle number 
(similarly, we will call this state particle cascade) can be 
found at this point by the formal substitution ui — 1, i.e. 
a = and 



d + n 



(27) 



Note that doing this substitution at later stages would be 
confusing since the explicit expression for /i also contains 
a. Note also that on turbulent states I[n] = 0, therefore, 
transport of all quantities except energy is zero for en- 
ergy cascade. For particle cascade, which describes Bose- 
condensation 83, 84], the transport of energy is zero. 

The reader should bear in mind that only those solu- 
tions that describe the transport of energy towards the 
ultra-violet, > 0, are relevant for the problem of ther- 
malization after inflation. The sign of fluxes for station- 
ary turbulence of three- and four-wave collision integrals 
was found in Ref. [s^. 

sign 5'' = sign [as{s — a)] . (28) 

In thermal equilibrium ncxw~^, i.e. s — a. Therefore, 
energy turbulence is directed towards the ultraviolet if 
the distribution function with increasing momenta falls 
off faster than in equilibrium, s > a. As we will see, 
in the X(f>* model this condition holds in d = 3, but is 
violated at d < 2. Therefore, we believe that simulations 
of the thermalization in this model at d < 3, see e.g. 
Rcfs. 52, 56, 58, 63], may not reflect all aspects of the 
physical problem of reheating after inflation correctly. 



C. Self-similar evolution 

In an analytical approach to non-stationary situations 
(e.g. when describing free turbulence) it is usually as- 
sumed that the evolution is self-similar js^ [s^ . As we 
have shown, the evolution is self-similar, indeed, at late 
times in our numerical integration of the 0^-model, see 
Section^ Below we consider self-similar substitutions in 



anticipation that they provide a valid leading description 
of thermalization in the class of models we consider. 

Let no{k) be a distribution function at some late mo- 
ment of time to, when the regime of self-similarity has 
been already established. The subsequent evolution can 
be described as rescaling of momenta accompanied by a 
suitable change of the overall normalization 



i{k,T) = A'^noikA) , 



(29) 



where we have defined t = t/to, 7 is some constant and 
A — A(t) is some time dependent function satisfying 
A{1) ~ 1. Both, A{t) and 7, are determined by the 
solution of the kinetic equation © . 

In some cases the collision integral may contain an ad- 
ditional explicit time dependence which can be isolated 
as an overall factor B{t). This factor may be induced by 
time-dependent classical backgrounds like the scale fac- 
tor of the expanding universe or the zero-mode of the 
inflaton field. It is convenient to rescale the collision in- 
tegral by some typical rate F, / = BTI, such that B and 
/ are dimensionless. We use B{1) = 1 as normalization. 

When Eq. holds, the factor A'' of each distribu- 
tion function, Eq. (|29|l . can simply be taken out of F and 
out of the collision integral, which becomes a functional 
of riQ. After that we can use Eq. 123|l with £, = A which 
gives 



Iik,T)^A-'("'-'^-^BrikA[no] 



(30) 



On the other hand, the l.h.s. of the kinetic equation 
can be written as 



n{k,T) = A'^-'^A (^-fno + C 



dno 



(31) 



where we have defined C, = kA. Using F as a separation 
constant, the kinetic equation can be split into two: one 
for the shape of the distribution function. 



and one for the dynamical evolution 

dr 



(32) 



(33) 



We will not be concerned with H32|l here and simply as- 
sume that it has some non-trivial solution. The general 
solution of l|33|) is of the form 



A^e-p 



where 



and 



e = ^ f B{T')dT' + 1 
p Ji 



1 



7(m — 2) — /i 



(34) 



(35) 



(36) 
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We fix scales using the condition Tto = p. For a time- 
independent background B, i.e. i? = 1, it than follows, 
that Q — T and eq. H34|) simplifies to 



A = T- 

We will discuss this case first. 



(37) 



Self-similar evolution in time-independent background 



Substituting ^ in ^ we obtain 

n{k,T) = T-'^P no{kT-P) , 



(38) 



In applications of turbulence theory to thermalization, 
this solution is most important. Let kc be the initial value 
of some characteristic momentum scale, e.g. the scale 
where most of the energy carried out by a self-similar 
distribution is concentrated. According to Eq. (j38|l , with 
time this scale evolves as 



k,{T) = fc,(l) 



(39) 



The exponent p determines the speed with which the 
distribution function moves over momentum space and 
therefore defines e.g. the time scale of thermalization. 
This is a reason why we will be interested mainly in the 
value of the exponent p. Eg. H36|l . In applications to ther- 
malization after preheating the energy is concentrated at 
low momenta initially and should propagate to high mo- 
menta. This means that solution Eq. H38|l is physically 
relevant for p > 0. 

The exponent 7, which enters Eq. H36(l can be fixed by 
specifying appropriate boundary conditions, which are 
specified below. 

a. Isolated systems If the wave energy is strictly 
conserved it follows that 



const = J d'^koJkn{k,T) 



This gives 



J = d + a . 



(40) 



(41) 



Similarly, for the evolution with particle number conser- 
vation one obtains j — d. Here we would like to stress 
the following subtlety. Clearly, a simple self-similar sub- 
stitution Eq. (|29|l cannot account for energy and particle 
number conservation simultaneously, while both quanti- 
ties are conserved in a number of systems. If this is the 
case, one should choose the integral which gives dominant 
restriction of n^, i.e. the energy for energy cascade (ther- 
malization) and particle number for the inverse cascade 
(Bose-condensation) . For the problem of thermalization 
of ultra-relativistic particles this gives 



1 



(d+l)(TO-2)-/i 



relativistic 
energy cascade 



(42) 



However, describing thermalization in the non- 
relativistic limit, ujk = m -\- k^ /2m with k"^ jm? ^ 1, we 
can neglect the kinetic energy with respect to the rest 
mass in the normalization condition (|4U|) . i.e. we should 
use 7 = d, as in the case of particle conservation 



1 



d{m — 2) — /i 



non — relativistic 
energy cascade 



(43) 



b. Driven turbulence. In our lattice integrations we 
have found that particle distributions as functions of k 
follow a power-law in the wake of a propagating energy 
front, rikir) = {b(T)/kY , with exponent s being in agree- 
ment with the theoretical predictions for stationary tur- 
bulence. Such behavior is expected [sBj for the regime of 
driven turbulence in the presence of a stationary source 
(and then &(t) — const). However, for the case of free 
turbulence we are not aware of any predictions. Here we 
consider consequences of such a behavior assuming gen- 
eral 6(t) (the case of constant b being a particular case). 

Considering distribution the function in the region of 
low momenta, nfc(r) — {b/ky — no{kA) we find 



b cx A''/'- 



-(1 — 



(44) 



i.e. the transport of energy through the inertial range is 
stationary if 
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(45) 



This generalizes the concept of stationary turbulence to 
a system without sink. (Notice, that this requires a sta- 
tionary source in the infra-red.) In this regime the total 
energy in particles has to grow linearly with time. Con- 
sidering the r.h.s. of relation H40|l with 7 = s we find 



Pt = l/{d + a - s) 



(46) 



where we denote the exponent p for the case of a station- 
ary transport as pt to distinguish it from the exponent 
which corresponds to an isolated system, pi. Substitut- 
ing explicitly the exponent s of the spectra of stationary 
turbulence, Eq. if^ . we find 



Pt = 



(m- 1) 



{d + a){m - 2) - ^ 



(m - l)pi 



(47) 



The latter relation could have been also found using 
Eq. ^ and Eqs. ^ with 7 = s. 

c. Non-stationary source. Let us consider the some- 
what more general situation and assume that the energy 
inputted into (or taken out from) the system of particles 
changes with time as E{t) = Eqt^ . Clearly, the isolated 
system corresponds to r = 0, while a stationary source 
corresponds to r = 1. We will now have 7 — {d + a) — r/p 
and 



1 -\- r{m - 2) 
{d-\-a){m-2) - ji ' 



(48) 
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2. Time-dependend background 

We now consider a time-dependent B in Eqs. (|34|l . 
(|35|l . As an illustration we choose B{t) — t^'^, which 
gives 



e 



1 



(1 



(49) 



Note, and this is important for the interpretation of our 
numerical results, that the linear approximation for small 
times, r ~ 1, gives ~ r, which brings us back to the 
situation considered in the previous subsection. 

The late time behavior, r 3> 1, depends on the sign 
of 1 — K > 0. If 1 — K > 0, the distribution propagates 
to the ultraviolet without bound, A{t) cx t^(^^'^^p and 



with 



In other words, at late times A 



p = (1 - k)p 



(50) 



for any boundary conditions discussed above in para- 
graph jlll C lal - IIII C 1 cl 

However, A{t) approaches a finite limit at r — > oo if 
1 - K < 



A{t = cxj) = 



1 



1 -p 



1 



(51) 



The propagation of particle distribution functions to- 
wards the ultraviolet is limited. This has important con- 
sequences for the thermalization of massive particles in 
the expanding Universe, as we shall discuss in more detail 
below. 

Expressions Eqs. and l|5n|l are the main re- 

sults of this section. They determine the speed of prop- 
agation of the particle distribution in momentum space 
for a specific models. 



IV. STATIONARY STATES AND SELF-SIMILAR 
EVOLUTION IN SPECIFIC MODELS 

Here we apply the general results of the previous sec- 
tion to a number of particular models of interest. 

First of all we have to determine the scaling exponent 
/i of d^l (see Eq. ()21|l '). The scaling of oj is different in rel- 
ativistic and non-relativistic regimes. This is accounted 
for differently in the argument of the energy conserva- 
tion (5- function (where in the non-relativistic regime a;(fc) 
is replaced by /2m) and in the 1/w factors of rela- 
tivistic integration measure (where oj is replaced by m). 
To make the discussion of relativistic and non-relativistic 
cases uniform, we move u out from the relativistic inte- 
gration measure and define the function C/(fc, qi) 



(52) 



In what follows we will assume that in a dynamically 
interesting range of wave numbers U follows a scaling 
law 



With this definition 



dVL{k,q,) = U{k,q,)5\k^,,q,^,) \[ 



d% 



and we find 



fi = d(m — 2) ~ a + f3 



(53) 



(54) 



(55) 



We calculate the exponents fj,, s, and p for two 
classes of models. The first one is characterized by k- 
independent matrix elements, the second one has no di- 
mensionful parameters. The scalar field models which we 
integrated on the lattice belong to the first class. In the 
absence of a zero mode in the relativistic limit in (3 -I- 1) 
dimensions they belong to the second class as well. 



A. Theory with fc-independent matrix elements 

For models with k-independent matrix elements the 
scaling of U is determind by the w's, and we have P = —m 
in the relativistic regime and P — in the non-relativistic 
case. Eq. gives 

/i = d{m — 2) — 1 — TO (relativistic) , (56) 
/i = d{m — 2) — 2 (non — relativistic) . (57) 

Substituting these expressions into Eqs. ll^ we find 
that in this class of models the exponents p do not de- 
pend on the number of dimensions. In particular, for the 
energy cascade in an isolated system we have 

Pi = l/(2m— 1) (relativistic) , (58) 
Pi = 1/2 (non-relativistic). (59) 

For TO = 3 and m = 4 Eq. H58|) gives p = 1/5 and p = 1/7 
respectively. 

Substituting Eqs. (EH, ^ into Eq. ^ we find the 
exponent s 



= d- 



m 



(relativistic) 



(60) 



TO — 1 

s = d (non — relativistic) . (61) 

In the non-relativistic regime both exponents, Pi and s 
do not depend on to. 



1. Three-particle interactions, relativistic regime 

Three-particle processes appear in the model when 
interactions with the zero-mode are important, see Ap- 
pendix [Bland Section irvTTTl 
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According to Eq. H58|) for m 

cascade propagates with 



Pi 



3 the front of the energy 



(62) 



regardless of of the number of spatial dimensions, d. For 
the exponent s of particle distributions in the inertial 
range in d = 3 we find 



(63) 



Both exponents coincide with what is observed in our 
numerical experiments. Note that the exponent s is ex- 
pected to appear in the case of driven turbulence. In 
the case of free turbulence the wake of the propagating 
turbulent front does not even have to be a power law. 
Nevertheless, we do observe a power law with the expo- 
nent s = 3/2 to a very good accuracy. This might be 
not a chance coincidence. However, in d < 3 the the- 
ory predicts s < 1, the spectrum falling-off with k more 
slowly compared to thermal equilibrium, and one can get 
a different shape of particle distributions in d < 3 (but 
we still expect the exponent p to be given by Eq. (|58|) '). 



B. Relativistic theory with dimensionless 
couplings. 

The model in d = 3 which we have simulated on 
the lattice belongs to the class of models considered in 
this paragraph. In d = 2 dimensionless couplings appear 
in the A(/)^ model. Dimensionless couplings are generic 
and this case is not restricted to scalar field models, there- 
fore we consider it separately. 

If the collision integral does not contain any dimension- 
full parameters, it has to scale with /i = 1 and we find 
for the exponent pi of energy conserving propagation in 
an isolated system, Eq. (|42|l 



1 



(d-f 1)(to- 2) - 1 



(64) 



For the physical case of d = 3 and for a 4-particle pro- 
cesses (which should dominate at late times in the models 
we have considered numerically, see below) we obtain 



Pi 



(65) 



Note that for d — 2 and m = 6 we have pi = 1/11, in 
agreement with Eq. H58() . For the exponent s of particle 
distribution functions in the energy cascade we find, see 
Eqs. if^ 



d + 2 



m - 



(66) 



C. Explicit time dependence in the collision 
integral 

The self-similar evolution is modified when an explicit 
time dependence is present. Below we consider two spe- 
cific models with explicit time dependence in the collision 
integrals which appear in the problem of reheating. The 
first one is directly related to the relativistic scalar model 
we have simulated on the lattice and time-dependence 
enters via the coupling to the zero-mode. The second 
describes thermalization of non-rclativistic particles and 
the time dependence is induced by the expansion of the 
Universe. 



1. Non-zero classical field 

Typically, oscillations of the inflaton zero-mode do not 
decay completely during the initial stage of parametric 
resonance. Moreover, if the resonance parameter is large, 
parametric decay stops early, when only a small part of 
the initial inflaton energy has been transferred to parti- 
cles . The remaining oscillating zero-mode serves as a 
source in our turbulent problem. This source acts via two 
different channels. The first one can be described as a di- 
rect decay into the resonance band(s). The other channel 
is m-particle scattering when one or more particles have 
zero momentum. These particles belong to the zero-mode 
(which is a Bose condensate). While the zero-mode and 
excitations with k ^ can be viewed as the same parti- 
cles but with different momentum, the formal description 
is different. The presence of the zero mode (po leads to 
new specific terms in the collision integral with reduced 
number of particles participating in the interaction pro- 
cess and different (and time-dependent) couplings. 

The simplest example is 2 by 2 scattering in the 
model when one of the incoming or outcoming parti- 
cles belongs to the condensate. These scattering pro- 
cesses can be modeled as an effective 3-particle interac- 
tion. The corresponding 3-particle collision integral can 
be obtained from the 4-particle one with the substitution 



Up 

UJr, 



(2^)3^(3) (p)02 



(67) 



This gives an explicit time dependence in front of the col- 
lision integral, B — (/)Q(r)/0§(l), and reduces the num- 
ber of integrations by one, m = 3. Alternatively, the 
3-particle collision integral in the background of a zero 
mode can be derived from first principles, see Appendix 

m 

The turbulent exponents for the 3-particle scattering 
without explicit time dependence (i.e. (f>o{T) = 1), are 
given by Eqs. and (|^ . Both agree with what is 

observed in our numerical experiments, see SecL^l We 
show in Sect. IVII that the collision integral in our lattice 
problem is dominated by 3-particle interactions. Ther- 
fore, Eq. (|63|l for the exponent s seems to be indeed 
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applicable for the system considered numerically. The 
question of applicability of Eq. (|62|l for the exponent p 
deserves special consideration because the amplitude of 
the zero-mode changes with time. 

During the initial stage, when the total energy in par- 
ticles is small compared to the energy stored in the zero- 
mode, we can consider the amplitude of oscillations to be 
constant and the source of turbulence to be stationary. 
However, distribution functions should then evolve with 
Pt = 2pi, see Eq. (|47|) . At late times on the other hand 
we cannot neglect the decay of the zero mode. Numer- 
ical integrations show that the amplitude of the zero- 
mode decreases as a power law, </>o(t) oc t~'^. At late 
times this gives p — + (1 — k)p, see Eq. (I5UII . Numeri- 
cally K = 2/3, however, the conclusion that p — 1/15 
would be incorrect. First, for completely decayed zero 
mode the 4-particle collision would dominate, leading to 
p = 1/7. Therefore, in our problem we should expect 
p > 1/7 ai all times. Second, the condition r ^ 1 is 
not fulfilled during our integration time. Indeed, we ob- 
served self-similarity for 3600 < rj < 10000, see Fig. O 
which corresponds to r < 3. For t w 1 the solution of 
Eq. lEHl for A{t) coincides with A = t^p, while at 

T = 3 it deviates by not more than 5%. Therefore, in this 
time interval A{t) sa r^^/^. Similarly, the quantity A'^ 
with 7 = 4 for 1 < r < 3 (energy conservation) is close 
numerically to t"'', where q ~ 3/5. Hence the indices 
of self-similar evolution obtained in Sec. ^ are explained 
by free turbulence driven by three-particle interactions 
in the background of zero-mode. 



2. Non-relativistic regime in expanding universe 



where k = 3/b for the 4-particle process in theory, 
i.e. K — 3 and k = 6 for radiation and matter dominated 
expansion respectively. This gives 



B{T')dT' = 



l-(6go7yo(r-l) + l) 

b{K - 1)HqT]q 



We see that in the limit t oo 



A{t = oo) 



1-f- 



&(k - l)Horio 



(70) 



(71) 



where p is given by Eq. (|43|) . The particle distributions 
cannot propagate to high momenta and are frozen out at 



fcc(r = oo) = 



fcc(l) 



fcc(l) 



A{t = oo) [b{K - l)HoVo]P 



(72) 



In the traditional discussion of thermalization of parti- 
cles in the expanding Universe, see e.g. 0, the expansion 
rate, Hq, is compared to the to the rate of interactions, 
which in our case can be identified with rjQ (see the nor- 
malization factor in Eq. (|34|l ) . It is concluded that parti- 
cles can not thernialize if Horjo > 1 while they can reach 
thermal equilibrium when i?o% < 1- Equation l|72|l tells 
us that thermalization is indeed impossible for i?o^o > 1 
since the distributions do not move towards high mo- 
menta in this case. However, it is not guaranteed that 
the equilibrium is reached even if Hor]o <^ 1. The sys- 
tem may thermalize only if fcc(T = oo) is not smaller 
than the typical values of momenta in eventual thermal 
equilibrium. 



Let us consider now non-relativistic particles in an ex- 
panding universe with physical dimension d — We will 
be working in the conformal reference frame, Eq. In 
these coordinates the expansion of the universe is simply 
accounted for by multiplying all bare mass parameters, 
M, by the scale factor. This is true both for the origi- 
nal field equations and for the kinetic equations (which 
are derived from the former). Factors of lu in the mea- 
sure Eq. should be replaced by Ma{r]). Therefore, 
in the non-relativistic regime the collision integral in the 
expanding universe can be obtained by multiplying it by 
the scale factor in some negative power. 

In conformal reference frame the solution of the Fried- 
mann equations for the scale factor as a function of 
T = Tj/riQ can be written as 

a'' = bHoVoir - 1) + 1 , (68) 

where Hq is the value of the Hubble parameter at time 
rjQ. For the radiation dominated expansion b = I, while 
6=1/2 for the matter dominated expansion. Hence, the 
function B{t) takes the form 



V. TWO INTERACTING SCALAR FIELDS. 
NUMERICAL RESULTS. 

In this Section we present the results of lattice cal- 
culations of reheating in the model of two interacting 
fields. As in the one field model presented in Section ITTI 
we again consider the massless case, for which the use 
of conformal transformation allows mapping of the dy- 
namics in expanding Friedmann universe into the case of 
Minkowski space-time. This permits a long integration 
time on a fixed lattice. 



A. The model 

At the end of inflation the universe is very close to a 
spatially flat Friedmann model. It is convenient to work 
in conformal coordinates where the metric takes the form 
ds^ — a{rjY {drf — dx^). We consider two scalar fields 
$ and X whose dynamics are determined by the action 
S = J dt d^x y/—g£{^, X) with Lagrangian density 



B{T)^[bHom{r~l) + ir'' ■ (69) 



^ = 25^"9^*a,$ + -g^^d^xd^x - y($, x) (73) 
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and potential 



-^X\ (74) 
4 



We identify the field $ with the inflaton. Therefore 
A$ ~ 10"-^^ Q, S I3- Inflation ends at time 770 when 
($(770)) =^0.35 Mpi. 

We use the following set of coordinate and field rescal- 
ings which bring the system into a dimensionless form 
suitable for numerical integration: 



dxo 



X 



dyi = dxi Ai/^A 



= $ A ^a{f]) 
X = XA-ia(r?) 



(75) 



(76) 



Re-scaling of the fields with 0(77) in Eq. (|76|) rotates the 
scale factor away and maps the model into a scalar field 
theory in Minkowski space-time. The classical equation 
of motion have two independent parameters 



5 = A$;ic/A$, h = Xx/M 



and simplify to 



(77) 



(78) 
(79) 



We choose A — {^{rjo)), so that the initial condition for 
the infiaton zero-mode reads {^{rjo)) = 1. The equations 
()78|l and (|79|l . however, are independent on the particular 
choice of A. At 77 = 770 all correlation functions of $ 
and X on subhorizon scales characterize a vacuum of 
fiuctuations around the inflaton mean value. 



B. Results of numerical integration 





FIG. 4: Particle distributions in the self-similar regime for 
h — lOg as functions of the corresponding wave kinetic en- 
ergies rescaled by the current zero-mode amplitude (j>o.. Up- 
per and lower panels correspond to x (j) fields respec- 
tively. In both cases from left to right the plots are taken at 
77 = 1000, 1500, 2000. 



We have studied the two-field model using the follow- 
ing set of couphng constants: A$ = 10^^"^, g = 30, 
and h was varied in the range O.lg < h < W^g. We 
will see below that different values of h lead to differ- 
ent duration and different relative importance of the spe- 
cific dynamical regimes, as it was already argued for in 
Sec. Ill B1 These are: the regime of parametric resonance, 
the regime of stationary (or driven) turbulence and the 
regime of free turbulence. These issues will be addressed 
later in this Section, which we start with the discussion 
of particle spectra. 

1. Spectra 

The particle spectra in the two field model at late times 
are very similar to what we have observed in the one field 



model and have the same turbulent exponents. Namely, 
in the inertial range nk is a power low with the exponent 
s ~ 3/2, for both fields x ^-nd (j>, see Fig. 0] And both 
fields evolve in a self-similar way with p = 1/5 at suffi- 
ciently late times, when the energy in particles became 
comparable to the energy in the zero-mode, see Fig. |3 
Both exponents, s and p, correspond to turbulence sup- 
ported by 3-particle interactions. 

There are some differences however. For the consid- 
ered range of parameters, the coupling of the excitations 
to the medium is rather strong, which induces large effec- 
tive particle masses, see Appendix IA3I Therefore parti- 
cles are non-relativistic already in the part of the inertial 
range. Namely, ~ 5.50 and ~ 1.70. This man- 
ifests itself as ~ fc~'^-power-law behavior, which is again 
consistent with domination of 3-particle interactions, see 



14 



xlO' 



3500 
3000 
2500 
2000 
'l500 
1000 
500 






FIG. 5: Spectral energy distributions for x (upper panel) and 
(j> (lower panel) in the model with h = lOg. In each panel 
we plot the wave energy per decade found in lattice integra- 
tions at three moments of time, rj = 1000, 1500 and 2000. 
In the lower-left corner of each panel are the same graphs 
transformed according to the relation inverse to Eq. 0- 



FIG. 6: Different regimes of the evolution of the x field for two 
values of self-coupling, h — lOg and h — lOOg. The dashed 
lines correspond to a linear growth of energy in the x field 
with time, oc rj. 



oscillations, ujq « M, tells us, see Appendix Second, 
the collision integral, Eq. ljB38p . being substituted into 
expression for the energy flux, Eq. 1)15(1 . will have ap- 
propriate universal scaling behavior in terms of kinetic 
energy, e^, but not in terms of k. Therefore, the kinetic 
energy is indeed the appropriate variable for the case of 
3-particle interactions in the presence of zero mode. 

For h > g the spectra look stationary in the inertial 
range after rescaling e by the current zero-mode ampli- 



tude (j}o ~ r]~^/^. This is similar to the one field case, 
(see fig. 0J. However, ioi h < g we found 0o ^ V^"^^^ ^ 
but the spectra still appear stationary after rescaling by 



This can be understood in the hght of Eq. 
6(t) = r^^/'^ is consistent with the choice 7 = 4, s = 3/2 
and p — 1/5. Hence, the decreasing amplitude of distri- 
bution functions in the region of low k simply reflect the 
energy conservation in the system. 



Eq. H59|l . This can be expressed as a single power law 
if particle distributions are plotted as functions of rela- 
tivistic kinetic energy, 

ek = ujk - M , (80) 

where M is the effective particle mass. Indeed, in the 
relativistic region we have Uk oc k^^'^ oc e^^^, while in 

the non- relativistic region we obtain n/j oc fc'^ oc e^^^. 
For this reason, the particle distributions were plotted in 
Fig. 01 as functions of e^. The particle distributions for 
the X field appear in this variable as featureless single 
power law. This can be easily understood. First, the 
energy transport for 3-particle interactions in the pres- 
ence of zero mode corresponds to the transport of kinetic 
energy, as energy conservation law in elementary scatter- 
ing process, which involves the frequency of zero-mode 



2. Stationary and free turbulence regimes 

Let us demonstrate now that the regime of station- 
ary turbulence does occur in the two field model. This 
regime is expected to appear in the case of large values of 
dimensionless parameters, g ^ 1, ft- ^ 1, when paramet- 
ric resonance stops early, while the total energy is still 
stored in the zero mode. 

We found that in the relevant range of parameters the 
description in terms of particles, which we were using so 
far, deteriorates. The reason is that in this language at 
large couplings there is no unique way to split the total 
energy density of the system into contributions coming 
from zero mode and fluctuation field. 

To deal with this problem we have quantified the en- 
ergy transfer in the following way. The quantity po = 
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FIG. 7: Spectral energy distributions at two moments of time, 
■q — 500 (dotted lines) and ri = 1400 (solid lines). We compare 
two models with different self-coupling, h = lOg and h = 
WOg. 

5 00 iVz ) gives a good measure of the total energy density 
stored in zero-mode oscillations, if it is measured at those 
moments of time, ry^, when the mean field crosses zero, 
4'o{ilz) = 0- III this way we can get rid of the ambiguity 
in accounting interaction energy between zero-mode and 
fluctuations. Similarly, we measure the energy density in 
the fluctuation field as = {x^)t and = {(j)'^)t for 
the X and (p fields respectively. Here (. . .)t means lattice 
and time averaging. We verified numerically that the 
sum of these quantities conserves with time and equals 
to the initial energy density. This is not true, however, 
when we measure the energy density in particles as LOkU}-. 
Both measures of particle energy converge at late times 
when the interaction energy becomes unimportant. 

This "kinetic" measure of the total energy density 
stored in particles as a function of time is shown in Fig.|Sl 
We compare models with two different values of h. Three 
different regimes are clearly seen in both cases. 

1. Parametric Resonance: The energy density p^ 
grows exponentially. This regime continues until 
re-scattering becomes important. The larger h is, 
the earlier resonance terminates. 

2. Stationary Turbulence: At later time the energy 
density in x particles grows linearly in time, which 
according to Eq. (O is a sign of stationary turbu- 
lence. During this period the energy density still 
stored in the zero-mode dominates the total energy 
balance. 

3. Free Turbulence: At some point the energy density 
in the zero mode drops below the energy density 
already stored in particles. Stationary turbulence 
cannot be sustained anymore and the regime of free 



turbulence, with conserved energy in particles, fol- 
lows. We may expect self-similar evolution of par- 
ticle distribution functions, which at late times are 
good quantities. 

In the model with larger self-coupling the parametric 
resonance stops earlier and only a negligible part of the 
inflaton energy is transferred to particles during the res- 
onance stage, see Fig El In this parameter range the 
transfer of energy from the inflaton into ^-field is dom- 
inated by a stationary turbulence. In the Sec. IVIIBl we 
show that if all coupling constants are of order of the 
inflaton self-coupling, the thermalization is a very long 
process and the Universe reheats to unacceptably low 
temperature, T ~ 100 eV. Therefore, some couplings in 
the sector of physical flelds (e.g. self-couplings, or cou- 
plings to the inflaton) in a realistic model have to exceed 
significantly the scale of the inflaton self-coupling. With 
larger couplings the thermalization proceeds faster. This 
is confirmed in our lattice integration, see Fig. [7| At 
earlier times the model with larger self-coupling contains 
less energy in x-particles, cf. curves ai rj — 500. How- 
ever, at later times this model takes over and the energy 
containing region moves faster towards ultraviolet in the 
model with larger self-coupling. 

With even larger self-coupling of the x-field, or its cou- 
pling to the inflaton, the period of stationary turbulence 
should become even more pronounced. In light of these 
findinffSLwe can also understand the results of earlier pa- 
pers IHQ. E.g., in Figures presented in Ref. 0, we 
see clear signs of driven turbulence, which was not iden- 
tified as such until now. In particular, in Ref. it was 
found that the energy in x-fluctuations grows with time 
as (X (Small deviation from cx t law can be due 

to the fact that the energy in zero mode decreases some- 
what and the source deviates from stationarity). This 
regime persists until the final integration time, when dis- 
tribution functions reach the boundary of the integration 
box, and even then the system is far from free turbulence 
regime. 

We conclude that in the models with an acceptable 
reheating temperature, the parametric resonance stops 
only when a negligible fraction of the inflaton energy has 
decayed. Therefore, in realistic models of the type con- 
sidered in the present paper, the major mechanism of 
energy transfer from the inflaton into particles is station- 
ary turbulence. 



VI. IS THE KINETIC APPROACH 
APPLICABLE ? 

In this section we confront the results of our lattice in- 
tegration with the predictions of kinetic theory and ad- 
dress the validity of the kinetic description at the ther- 
malization stage during our integration time interval. 

The particle distributions in the inertial range, n{k) 
k^^ with s K, 3/2, which we observe in the lattice 
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simulations, can be understood as corresponding to the 
scale invariant energy flux for 3-particle interactions, see 
Eq. H63(l . The observed exponent p — 1/5 of the self- 
similar scaling of free turbulence, can also be in accord 
with 3-particle interactions, see Eq. H62|l . However, in 
our case bare 3-particle couplings are absent and appear 
effectively in interactions with zero-mode. Therefore, the 
3-particle collision integral is multiplied by the amplitude 
of zero mode squared. Since the amplitude of the zero 
mode oscillations decays, one can expect p = 1/5 only 
during a small time interval, see Sect. IIII C 21 

Can 4-particle interactions be responsible for the ob- 
served scalings ? For 4-particle interactions pi = 1/7, 
see Eq. (|^ . which is not that far away from the lattice 
results, especially if one takes into account energy influx 
from the zero-mode. However, for particle distributions 
in the inertial range one should expect s — 5/3, which 
is not in a good agreement with the observed value of 
s = 3/2. Further, in view of Eq. H67|l one should expect 
the dominance of 4-particle scattering during the time in- 
terval when the variances of fluctuations are larger than 
(pQ. This is not the case during the time interval encom- 
passed by the lattice simulations, see Fig. ^ 

The outlined difficulties may give an indication that 
the weak turbulence description is not applicable in our 
case. In view of the importance of the issue, we per- 
formed a detailed study of collision integrals, anomalous 
and higher order correlators, as measurements on the lat- 
tice, and compared these with predictions and assump- 
tions of kinetic theory. 



A. Collision integrals 

To verify the extent of agreement between kinetic the- 
ory and lattice calculations, and to find out which pro- 
cesses dominate the collision integral in our problem, we 
carry out the following procedure. First, we numerically 
calculate the collision integrals using standard expres- 
sions, Eqs. lf7|l- HlU|) . and the particle distribution func- 
tions nk{r]) extracted from our lattice calculations. Sec- 
ond, using lattice data we calculate time derivatives of 
the distribution fimctions to see if the relation rik = Ik [n] 
holds. We limit ourselves to 3- and 4-particle collisions. 

The general relations, Eqs. (|7|)- (|10|l . for 4- and 3- 
particle collision integrals can be reduced to two and one 
dimensional integrations respectively, if the distribution 
functions are isotropic. Explicit expressions are given in 
Appendix m, Eqs. l|538l) and l|539)l . 

The numerically calculated values of J^'^-' and ij^^^ colli- 
sion integrals are shown in Fig.|Slin comparison with hk- 
Note that the collision integrals and hk take positive and 
negative values. For clarity we show only absolute values 
of these functions and indicate schematically the bound- 
ary between regions where hk is negative and positive. 
Roughly, in the inertial range hk is negative (recall that 
in this region the particle distributions can be approxi- 
mated as nk{ii) = {(po/ky and are decreasing functions 
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FIG. 8: Absolute values of h{k) and of ij^^ and /^^' collision 
integrals at r; = 5000. To the left of the arrow h{k) and colli- 
sion integrals are negative, while to the right they are positive. 

(3) 

Occupation numbers, n^, are also shown for comparison. 
agrees with h{k) to the left of the vertical dashed line. 



of time), while hk should be positive at larger k where 
the cut-off starts (recall that energy is flowing into this 
region). 

(3) 

We find that gives a reasonable approximation to 
hk practically in all range of k which is dynamically im- 
portant, which is to the left of the vertical dashed line in 
Fig.|Sl One reason for the disagreement between h{k) and 

(3) 

If. at larger k could be due to the fact that on the lat- 
tice some of the allowed resonant wave interactions of the 
continuum limit are not present ( cf. |87|). In any case, 

(3) 

in the region where /j, and hk disagree, the occupation 
numbers are relatively small, Uk < 10^, and this region 
should not contribute to the dynamics significantly. 

The ij^^ collision integral is about an order of magni- 

(3) . 

tude smaller compared with If, and is subdominant in 
the evolution of Uk, except on the very tail of the distri- 
bution, see Fig. |H1 The agreement between hk and /^^^ 
in the region of the tail is not coincidental - we observe 
it at all rj. 



B. Anomalous correlators 

Usually, kinetic equations are derived under the as- 
sumption (akCLq) <C {a*f,aq). However, this condition not 
always holds. For example, in the case of particle cre- 
ation by a time-varying classical background (e.g. in the 
region of parametric resonance) 

nfc = — Re(afe), (81) 
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However, in many cases it is necessary to trace the 
events further in time, e.g. to find out when and how the 
symmetry which was restored during preheating gets bro- 
ken later on, or to trace which fraction of baryon or Icpton 
asymmetry survives in the process of thermalization. Fi- 
nally, one needs to know when thermal equilibrium will 
be established. This gives e.g. the abundances of partic- 
ular dark matter particle candidates and other, possibly 
cosmologically "dangerous" relics like the gravitino. 

The explicit time dependence of the particle distri- 
bution functions, and the knowledge that the evolution 
is self-similar, n(fc,r) — nQ(kA) = t~'^p no{kT~P), 
which we have found in the present paper, may be use- 
ful here. Below we discuss some applications, limiting 
ourselves to field variances and to the problem of ther- 
malization. 



FIG. 9: Occupation numbers and absolute values of at 
ri = 5000 (dashed lines) and = 10000 (solid lines). 

where 

(flfeaq) = CTfc (5(k + q) , (82) 

see Appendix IbI In this case the anomalous correlators, 
tjfc, can not be neglected, since ak ~ Uk- This holds in 
general: if coherent processes are important, the correla- 
tors Eq. (|82|l may modify the dynamics of Uk- If this is 
the case, they should be included into the kinetic equa- 
tion. Since ak were neglected in the kinetic equations, 
Eqs. (|7|)- (|10|l . it is important to verify if the condition 
|cfc| <C rifc holds in our simulations. 

The correlators ak are shown for several moments of 
time in Fig. |51 In the inertial range the anomalous cor- 
relators are small indeed, \ak\/nk « 3 • 10~^, while this 
ratio is an order of magnitude larger in the region of the 
resonance peak {k « 0.5 at late times), which is expected 
behavior. The ratio \ak\/nk is growing also in the region 
of large k, reaching the value of 0.1 at A; = 8 at late 
times, see Fig. |31 To avoid confusion, note that k = 8 
corresponds to fc/0o ~ 25, which is the variable used in 
Fig.|5| We do not know if the growth of \ak\/nk at large 
fc is a lattice effect, but we can conclude that the kinetic 
equations in its simple form, Eqs. l(7|l- H10|) . should be 
applicable in the inertial range. 

VII. PHYSICAL APPLICATIONS 

Many different effects may occur during the stage of 
preheating. Some of these were discussed in the Intro- 
duction section. They have a common physical origin: 
rapid particle creation and large accompanying fluctua- 
tions of the classical fields involved. These findings are 
unaffected by our results, even in the case when only a 
relatively small fraction of the inflaton energy is trans- 
ferred to fluctuations during the initial stage of paramet- 
ric resonance. 



A. Field variances 

In some applications, basic observables like field vari- 
ances may already give the answer to the problem in 
question. This applies to the problem of symmetry 
restoration. To illustrate this, let us consider the Higgs 
field which is coupled to a x-field. In the vacuum state 
without condensate the mass squared of the Higgs field 
would be negative, — and the corresponding symme- 
try is broken. In the presence of the background of x~ 
particles, the mass gets "dressed", = — /x^ -I- gix^)- 
If the field variances are sufficiently large, the symmetry 
is restored (and is broken when (x^) < /i^/i?). 

If anomalous correlators are negligible, the field vari- 
ances can be calculated using expression 

varix) ^ {X')~{X? = / 7^^ - • (83) 

J (27r)'' LOk 

With the help of the self-similar substitution, Eq. H29|l. 
we find 

mr(x,T) = AT-'^+"mro(x) . (84) 

Here, the left hand side is taken at conformal time rj, 
while waro(x) on the right hand side is the variance at 
some earlier time rjQ. 

1. Relativistic regime 

a. Free turbulence. In this case 7 = d + a, see 
Eq. H41(l . and we find with a = 1 

var{x, t) = A^" varo{x) = r'^^' varo{x) ■ (85) 

For systems that we have studied numerically, pi — 1/5 
at early times which span the integration period. There- 
fore, in the free turbulence regime, we should expect 
var{x, t) = T^^/^ z;aro(x). This is in agreement with the 
results of our numerical integration, see Fig. ^ For late- 
time evolution, when 4-particle interactions will start to 
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FIG. 10: Time dependence of the variance of x field in the 
model with h — lOg considered in Sec. [V] 



dominate, we predict a slower decrease of the variance, 
var{x,T) = T"'^/'^ varo{x). 

Note that these resuhs have to be divided by the scale 
factor squared if expansion of the Universe is important. 

b. Driven turbulence. In the regime of stationary 
turbulence without sink, according to Eq. H45|l . 7 = s, 
which gives var{x, t) = ^s-d+a y(irg(^^y Using Eq. (|6U|I . 
we find s — d + a = l/(l — m) and 



var{x,T) 



_ -J)t/(l-m) 



'uaro(x) , (86) 



where relation Eq. (|47|l was also used. Therefore, 
during driven turbulence the variance should grow as 
war(x, r) = r^/^ wa?'o(x) in the models which we have in- 
tegrated numerically. This is indeed the case as Fig. [TUI 
shows. The transitional period from the regime of para- 
metric resonance to the regime of stationary turbulence 
at 77 ~ 10^ is slightly more pronounced in Fig. ^|as com- 
pared to Fig. This may be explained by the fact that 
different regions of momentum space are emphasized in 
and in var{x)- 



2. Non-relativistic regime 

In the case of free turbulence we have "f — d, and in 
Eq. (|84|l we have to substitute a = 0, which corresponds 
to LJk M in Eqs. lO, Therefore, var{x,T) = 

const. For driven turbulence we have ^ — s = d, see 
Eq. H61() . and Eq. (|84|l again gives var{x,T) = const. 

We see that in the regime of driven turbulence vari- 
ances are slowly changing functions of time (oc r^/^ in 
the relativistic case and war(x, r) = const in the non- 
relativistic case), while energy in particles grows fast, 
cx r in this regime. This is in accord with the fact that 
variances can be large right after the initial parametric 



resonance stage, while the amount of energy transfered 
during this stage is low and all energy transfers occur in 
the regime of driven turbulence. 



B. Thermalization in the absence of zero-mode 

We now apply the results obtained earlier in this paper 
to the general problem of thermalization of relativistic 
and non-relativistic scalar particles, both in Minkowski 
space-time and in expanding Friedmann universe. We do 
not restrict ourselves to the models which were studied 
numerically. Our analysis will be based on expression 
(|29ll with the factor A{t) being specified for a partic- 
ular modeled. This expression describes a self-similar 
propagation of the distribution functions into the ultra- 
violet. In a classical theory this evolution continues with- 
out bound (unless we consider a non-relativistic theory 
in expanding universe). 

The classical evolution stops when a system reaches the 
quantum regime where it can relax to the Bose-Einstein 
distribution. We adopt that this happens when in a re- 
gion of momenta, kf, which saturate the energy integral, 
the occupation numbers became of order one. In this 
subsection we consider the case of free turbulence. Then, 
one can estimate fc/ using energy conservation and ap- 



proximating the energy density as p 



in the 



region where ^ 1- On the other hand, initially the 
energy was deposited into particles with lower momenta, 
denoted below as ki. The relation hi ~ M^, where 
is the inflaton mass, determines the scale of initial mo- 
menta. Eq. H39|l gives for the time needed to thermalize 
a system: 



„th 



{kf/k, 



(87) 



Actually this should be considered as a lower limit on the 
thermalization time since we have to add a time which 
the system will spend in the quantum regime. 

As an idealization of the thermalization process we 
consider the evolution of a sub-system of excitations of 
a field Xj assuming that a fixed part of energy was de- 
posited into it initially, while since then x evolves as an 
isolated system. In this subsection, for estimates of the 
thermalization time we neglect the presiding regimes of 
parametric resonance and of stationary turbulence, since 
they are much shorter if the relevant coupling constants 
are not drastically different. We consider the possibility 
of (partial) thermalization in the regime of driven turbu- 
lence in the following subsection. 

As a first step we will find the thermalization time 
which follows from the exact self-similar solutions ob- 
tained above. Then we will show that in all cases which 
we consider, the result coincide, parametrically, with the 
"naive" perturbative estimates. Doing this comparison 
we neglect all numerical coefficients. 
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Relativistic regime 



Equation (|57|l gives 



The expansion of the Universe is easily treated in confor- 
mal reference frame. We have pf = pi = c^ptotj where c-^ 
is the fraction of the inflaton energy deposited into the 
field X during preheating and driven turbulence. This 
finalizes the answer. The result is general and is valid 
for any model. The initial inflaton energy can be written 
as ptot ~ ^'i4>t ~ where 00 is the initial am- 

plitude of inflaton oscillations. We find with p = 1/7, 
Eq. (|(j5|l . which corresponds to a relativistic theory with 
dimensionless couplings 



_th 



^7/4 



Mpi 



7/2 
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(89) 



We used here the inflaton parameters, Af^ « 10 ^Mpi 
in M^cjP model, or = \/\Mpi with A « 10"^^ in the 
A(/)^ inflationary model. 

To avoid confusion, note that the definition of r is 
different in different models since it involves to ~ F^^. 

We can assume c^,- « 1 if the Universe expansion can be 
neglected (this may be of interest for problems outside of 
inffationary cosmology) , and if the number of competing 
channels (other fields beyond x to which the initial energy 
can be deposited) is not large. 

a. Minkowski space-time Let us show that Eq. I|88fl 
agrees with the "naive" perturbative estimate. For this 
estimate we define r as r = F, where tn is the per- 
turbative estimate of the thermalization time t~^^ ~ an. 
In the Xcj}'^ model a - X'^/T^, n ^ T^, T ^ kj and 
therefore t~^^ ~ ^'^^f- On the other hand, parametric 
resonance stops when the rate of re-scattcring from the 
resonance band becomes equal to the rate of particle pro- 
duction fj, ~ ~ ki. This gives F ~ fc^ and we find 
T ^ rtu ~ ki/X^kf. Now, p ^ kj and p ^ kfnk, where 
Uk correspond to the typical occupation numbers at the 
time when parametric resonance stops, n^ ~ 1/A. We 

obtain tfjT ~ (^p^^'^/ki) , in agreement with Eq. (|5S|) . 

b. Friedmann universe In this case we can estimate 
the final temperature as T ~ kf/a^r). Let us consider 
a radiation dominated universe with a(T) = iJo77o(T — 
1) -t- 1, see Eq. We neglect the rapid epoch of 
stationary turbulence, and 770 corresponds to a time when 
the evolution of x is driven by its self-interaction with self 
coupling A^, i.e. r/^^ - F - A^n^fc, - X\{c^p^^^f / kj , 
where we have used p^ ~ kfnk- On the other hand ~ 
p[ll /Mp\. Combining this with Eq. H88II we find a(r) — 



Hor]oT 



-1/4 1/4, 



Ptot /'^iMpi. For the final thermalization 



temperature we obtain T ^ kf /a{T) 



we have used kf ^ p^ 
estimate, an ~ H. 



1/4 



Cx^X^AIpi, where 



This again agrees with the naive 



NumericaUy, T - A^Mpi ~ 100 eV, if we use the 
strength of the inffaton self-coupling, A « 10""'^^. Ther- 
fore, in a realistic model, at least some couplings should 
be significantly larger than this scale. 



2. Non-relativistic regime 

Now we consider the thermalization of X-particles of 
mass Mx in the non-relativistic regime. We assume 
that the relaxation is due to the self- interaction XxX^. 
The particle number conserves in the conformal reference 
frame in this regime, and Eq. (|87|l gives 



^th 



,1/3P 



1 



Pf 

Mx 



1/3' 



1/p 



(90) 



where cx is the fraction of the infiaton energy which ini- 
tially was deposited into the field X (this fraction should 
be measured at the time when the self-similar evolution 
starts). In the present case p — 1/2, see Eq. (|59|l . and 
similarly to Eq. I|89(l we find for the relaxation time 



Cx 



M^Mx 



2/3 



Cx 



M 



-1 2/3 



X 



(91) 



In this expression AI^/Mx ~ 1 since Bosons which are 
much heavier than the inflaton are not created, and in the 
opposite regime X-particles would have been relativistic. 

As we have seen in Sec. IIV C 21 there is no real relax- 
ation of massive particles when the expansion has become 
important. If some relaxation happens, it should occur 
during the time interval when the scale factor does not 
deviate significantly from its initial value. Then the ex- 
pansion can be neglected and the relaxation proceeds as 
in Minkowski space-time. Let us show that the expres- 
sions above agree with the "naive" perturbative estimate 
in the latter case. 

a. Minkowski space-time The perturbative relax- 
ation time in the final state can be estimated as t 



t 



/M^ and n ^ 
On the other hand. 



jL3 



R 



Therefore 
the rate in the 



van, where a ^ 
-^^^^Xlk^Ml 
initial state is given by a similar expression, but is mul- 
tiplied by large occupation numbers in the initial state 
|88| (which can be viewed as Bose-amplification factor), 
F ~ vannk ^ X\n'^/kfM\, where we used n ^ kfn^- 
We obtain t^F - n'^/kjkf - {p/MxY'^/kf, where 
p — Mxn. This agrees with Eq. (|^ . 

b. Friedmann universe To estimate the thermaliza- 
tion time and temperature we need to know the typical 
rate of reactions and the value of the Hubble parameter 
at the beginning of self-similar evolution. For definiteness 
we consider the situation which arises after preheating in 
the massive inflaton model coupled to a heavy fleld X. 
We assume that self-coupling of the X-fleld is sufficiently 
large, such that the "parametric" decay of the inflaton 
is halted by X-rescattering on each other. Using results 
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of Ref. [8j we conclude that at the moment when the 
inflaton zero-mode decays completely, the energy den- 
sity in the X-field can be estimated as ~ Af^/A^, 
while the rate of re-scattering is r^^j"^ ~ Mx- This gives 

fl"o??o ^ Xx^^'^Mx/Mpi ^ (?~^^^, where q is the initial res- 
onance parameter. Since q can be very large, the product 
HqI]^ can be small and the expansion is not significant 
at the initial stage of the self-similar evolution. On the 
other hand, the time needed to reach the quantum regime 
is of order r'^ {kf /kiY^P . Since particle number con- 
serves during the period of self-similar evolution we have 
~ M^/Xx, while ki ~ Mx at the end of parametric 

resonance stage. This gives r*'* ^ (1/A^)^/'^p ^ X^^^^. 

The condition HorjoT*'^ < 1 gives X^^ > Mx/Mpi as a 
necessary condition to reach a thermal state before the 
freeze-out of distribution functions. Using inflationary 
normalization, we conclude that non-relativistic particles 
created in "parametric resonance" have a chance to ther- 
malize between themselves in an expanding universe if 



On the other hand, the energy in the subsystem of 
X-particlcs grows in the regime of driven turbulence as 
Px{t) = T Px{l), and should not exceed the total energy 
stored in the inflaton zero-mode oscillations. The initial 
energy can be estimated as Pxi^) ^ ^ros'^x' where k^cs ^ 
q^^^uj^ and q is the resonance parameter: q = A0^$q/M| 
in the M^$^ inflaton model, or g = A^^/A^ in the A^'' 

inflaton model. This gives Px{^)/ptot ~ ^(pxl^x^ ^^"^ 
obtain the bound 

r < Xx/X^x ■ (93) 

We conclude that the quantum domain can be reached 
in the regime of driven turbulence if A^ > A^^^'*''^^'' — 

^i/x ^ lO^''- Here we have used s = 3/2, p — 2/5 
and A^^ ~ lOA^. These values are realistic, therefore, 
physical implications of driven turbulence in applications 
to thermalization deserve further study. 

VIII. CONCLUSIONS 



C. A faster route to thermalization ? 

Considering the regime of free turbulence, we have ob- 
tained estimates for the thermalization time which are 
in agreement with "naive" perturbation theory. It was 
important in these estimates that the relativistic free 
turbulence propagates with p = 1/7. This should be 
true at sufficiently late times, when all effects related to 
zero mode become insignificant. However, free turbu- 
lence driven by 3-particle interactions in the presence of 
a zero mode evolves with p — 1/5. The evolution of the 
front of particle distributions is even faster in the case 
of driven turbulence, when p = 2/5. If the quantum do- 
main is already reached during one of these two stages 
our estimates for thermalization should be changed. 

Here we consider the question whether a subsystem of 
X-particles can reach the quantum region in the regime 
of a stationary turbulence. 



1. Driven turbulence 

The quantum domain is reached in the regime of driven 
turbulence if the power law of the inertial range will ex- 
tend up to rifc 1. In other words, rife = {k/kT)~'^ should 
be valid up to fc = fcy- Let us consider the model were the 
largest coupling is the self-coupling of the x-field. The 
normalization of Uk can be fixed if we recall that in the 
region of the source, k ki, the X"P9'rticle distribution 

is given by Ux 1/^x- This gives kr ~ kiX^^^^ , or the 
time needed to reach the quantum region is given by 

r ^ A-i/^f , (92) 

where we have used Eq. (|87|) . 



We have studied the process of thermalization of clas- 
sical systems, which at some point in their evolution are 
in a highly non-equilibrium state with energy being con- 
centrated in a deep "infra-red" region of momenta. Such 
states naturally appear e.g. during reheating of the Uni- 
verse after cosmological inflation. We have shown that 
the process of relaxation in such systems can be divided, 
in the general case, into three distinct stages. 

In the models of the type we have considered in this 
paper, the initial stage of preheating is powered by 
parametric resonance. During this initial linear stage the 
rate of energy transfer is the fastest. The energy in par- 
ticles grows exponentially. However, in the physical situ- 
ation of reheating after inflation, the coupling constants 
have to be sufficiently large to insure an acceptably short 
time-scale of the subsequent thermalization, while with 
large couplings, only a negligible fraction of the initial 
inflaton energy is transfered into fluctuations during the 
parametric resonance stage H, @ . 

We have shown that in such situations the linear stage 
is followed by the regime of a driven stationary turbu- 
lence. During this stage, the energy in particles grows 
linearly in time. The regime of stationary turbulence 
stops as soon as the energy in particles starts to domi- 
nate the overall energy balance. Therefore, this regime 
is a major mechanism of energy transfer from the oscil- 
lating inflaton zero-mode into other species in realistic 
models of the type we have considered here. This period 
of evolution is also prompt. It should be noted that the 
source which drives the turbulence is powerful because 
coherence effects are still strong in the relevant region of 
momenta. 

The subsequent long stage of thermalization classifies 
as free turbulence. This stage should be generic. The 
energy in particles is conserved during this epoch, while 
the shape of the particle distribution function changes 
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in a self-similar way with the front of the distribution 
propagating into the ultra-violet. This stage continues 
until the quantum regime is reached and particles can 
relax to Bose-Einstein distributions. Applying conven- 
tional kinetic theory we have calculated analytically the 
time needed to equilibrate a system and the resulting 
temperature in terms of coupling constants and initial 
inflaton amplitude. The result coincides parametrically 
with the "naive" perturbative estimates |3] . 

We made a comparison of kinetic theory with the nu- 
merical integration of scalar field models on the lattice. 
We show that, at late times, the kinetic approach is ap- 
plicable, resulting in a weak wave turbulence regime pq . 
In the models considered numerically, the evolution is 
driven by three-particle scattering in the background of 
zero-mode oscillations. The characteristic exponents cal- 
culated within the framework of wave kinetic theory are 
consistent with the results of our lattice simulations. 
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APPENDIX A: NUMERICAL PROCEDURE 

In our study we have developed and employed a higher 
accuracy version of the LATTICEEASY code 89]. Var- 
ious correlation functions were measured with the use of 
Fast Fourier Transform (FFT), adopted from Numerical 
Recepies Essential details of our procedure can be 
found in Refs. 0, 0) • Here we describe specific 
choices of parameters, significant important differences 
in the integration scheme and give exact definitions of 
lattice observables. 



The vector index e runs over the three orthonormal direc- 
tions of the lattice. Al is a fourth order approximation, 
i. e. (A — Ai^)^{x) ~ 0{b'^) for a differentiablc function 
^{x) . The Fourier transform of Al differs from that 
of A, which would be given by multiplication with k^. 
Therefore the dispersion relation for a massless field on 



The numerical integration was done on a 3-D cubic 
lattice with periodic boundary conditions. The lattice 
is parameterized by the box size L, and the number of 
lattice points per dimension, N. These give the lattice 
spacing b = L/N and the total number of lattice sites, 
in three dimensions. 

The results presented in the paper are taken from sim- 
ulations with 256"^ lattice sites and a box size L chosen 
to fit a particular problem. For example, in the case of 
Eq. 0, i = T.Stt. With this box size the infrared modes 
which belong to the resonance band are still well repre- 
sented, while the ultraviolet lattice cut-off is sufficiently 
far away from the occupied modes, therefore the particle 
spectra are not distorted even at late times. We have 
studied the dependence of our results on the lattice- and 
the box size to avoid lattice artifacts. 

The finite-differences scheme that was used is 2nd or- 
der in time and 4-th order in space. 

a. Finite- differences scheme 

We write the equations of motion Q or (|78|l . (|79|l as 
fourth order finite differences on a three-dimensional spa- 
tial cubic lattice with periodic boundary conditions. The 
corresponding equations were evolved with the use of a 
symplectic integration scheme. Details are as follows. 

Particle wave numbers are discrete on the lattice, 
k — (rii, n2, n3)fco: where —N/2 < nj < N/2 and fco = 
{2t:)/L. The phase space is restricted to fco < fc < fcmax, 
where fcmax = V^koN/2. To avoid distortion at high mo- 
menta, it is desirable to take large N. This, however, is 
limited by the capabilities of the computer used. The 
choice of small values for L is also prohibited since that 
will lead to infrared distortions and may even move the 
resonance band out of the integration box. The problem 
is alleviated by the choice of a finite-differences scheme 
which is fourth order in space. This can be quantified in 
the following way. 

The lattice realization of the Laplacian in our scheme 
is given by: 



(Al) 

I 

the lattice is also different, and is given by 

'^L(fc) = ^Y.(l^l '^o^^^') + I cos(2afc,)') , (A2) 

i— 1 ^ ^ 

We find that < fc^ and ul - 0{k%'^) for smaU 
fc. Numerically, for fc < fc,nax/3, the relative difference 
between fc and is less then a percent, while for larger 
fc it grows up to about 30 percent difference at fc = fcmax- 



1 4 5 

— <I>(a; + 2be) + -<^{x + be) - -<^{x) 



4 1 
-$(a; - be) - -^'^{x - 26e) 
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FIG. 11: Deviation of the dispersion law for massless excita- 
tions on the lattice, (fc — LJL)/k, for second and fourth order 
finite-difference schemes. 



This means, that we can expect essentially undistorted 
self-similar and turbulent solutions on the lattice, if the 
dominating modes have wave vectors k < fcmax/3. In 
the case of a second order realization of Al, we find a 
considerably smaller available phase space, k < fcmax/lO. 
This is illustrated in Fig.llllwhere we plot {k—LOL)/k as a 
function of k for the second and fourth order calculation 
schemes on the lattice L = 7.57r, used in our integration 
of the problem Eq. Q. We see that up to k = 10, 
which essentially encompasses the support region of the 
distribution functions, see Fig. |31 the dispersion law on 
the lattice represents the continuum correctly. That is 
why self-similarity was not distorted on our lattice and 
could have being detected. (The small deviations from 
self-similarity, which can be observed at the very tail of 
the distribution and at the latest time, see Fig. 13 are 
caused by the distortion of the dispersion law which starts 
to be non-negligible here). 



Classical Approximation and Stochastic Initial 
Conditions 



The initial linear stage of parametric resonance has a 
complete quantum description, which is best expressed in 
the language of Bogoljubov transformations. The quan- 
tum description of this linear problem can be mapped 
into an equivalent classical problem |^. In our dimen- 
sionless variables the initial conditions for the classical 
description are given by the following probability distri- 
bution for field fluctuations in Fourier space: 



P[ip, ^] ~ exp 



d'k2Lot{To)\i^k 



(A3) 



Here "■0" should be replaced by one of the fields </) or x 
that are the dynamical variables in the simlated equa- 
tions lig or |(^), (O, and (ryo) 
The symbol 5^ ( 



9- 



\Jk\ -)- 3, while 
) is a functional 



Dirac-distribution so that the canonical momenta are 
locked to the classical trajectory. 



c. Measured Quantities: 

We measure various physical quantities both in con- 
figuration space and in Fourier space. In configura- 
tion space it is convenient to measure the zero mode, 
(f)Q = {(/)), and the variance, var{if)) = {(j)'^) — 0g. In 
Fourier space we measure particle number and other cor- 
relators. 

For large N lattice averages basically coincide with the 
statistical ones (ergodic theorem). We use this fact to 
measure expectation values (zero- modes), variances and 
higher cumulants of fields and their conjugate momenta. 

Spatial Lattice Averages: For averages defined in con- 
figuration space (O) = J d?x O, which on the lattice 
is expressed as the sum over the lattice points (O) = 

Fourier Spectra: For monitoring purposes we make a 
FFT transform at least every period of inflaton oscilla- 
tion. The wave amplitudes of fourier transformed fields 
are defined by Eq. ljB20|l . see Appendix IbI In the dimen- 
sionless units that we use in the numerical simulation the 
physical wave amplitudes take the form 



H = 



1 uj'lil)k + iipk 



(27r)3/2y^ 



(A4) 



where again "-0" stands for the dynamical variables 4> 
or X in the equations Q or (|75|l . ifT^ . The dimension- 



where 



Mak- 



less frequencies are given by = \l k\+ rri^g 

ing use of ak, we calculate various correlators, n(fc) = 
(aj^afc), cr(fc) = (afca_fc), (a*a*aa), etc. The first one, 
which corresponds to the particle occupation numbers, is 
of prime interest. 

Note that with this simple definition of quasi-particles 
the Hamiltonian is not diagonal in terms of and 
af. wave amplitudes if interaction energy is important. 
Therefore, the related definition of e.g. particle number 
is good only for modes with dominating kinetic energy. 



APPENDIX B: KINETIC EQUATION FOR 
CLASSICAL WAVES 



Following the general approach of Refs. p&L 1801 we 
derive the wave kinetic equation for the classical system 
of interest, the massive Xif)*- theory in d dimensions with 
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Hamiltonian density 



i(V0)^ + ^<^2 + ^0S (Bl) 



and in the presence of an oscillating classical background. 
Wc assume random wave fields which are statistically uni- 
form, i.e. the equal time correlation functions of cj) and 
its canonical conjugate momentum (j) are homogeneous 
and isotropic. We also assume the field to be weakly 
interacting. 

The first step in the derivation of the kinetic equation 
for an arbitrary system is to find Fourier wave ampli- 
tudes, afc, such that the quadratic part of the Hamilto- 
nian is diagonal in afc, i.e. : 



Ho 



(B2) 



The general equation of motion for the wave amplitudes 
is 



da{k, t) 
Jt 



da{k,t) . . SH^nt 

lUJkayk, t) — i- 



dt 



Sa*{k,t) 



where Hint = H — H2. The first term in the l.h.s. is due 
to a possible explicit time-dependence in the definition 
of Ofc, which can appear for example in case of a time- 
varying background. 

In the kinetic approach we want to get rid of rapidly 
varying phases of the wave amplitudes, i.e. to derive the 
equation for the slowly changing "occupation numbers" , 



~ a^a^.. To achieve this we multiply Eq. (|B3|I by a^, 
subtract the complex-conjugate expression and average. 
The result will contain higher order correlators induced 
by interaction terms. The resulting BBGKY-hierarchy of 
equations for Fourier cumulants can be solved, e.g., in the 
random phase approximation in consistent perturbative 
expansion. 

In the case of (jBip the wave amplitudes for the fluctu- 
ation fields are solutions of 



(2^ 
A/2wfe 
(2^)''/2 



V2i 



(B4) 
(B5) 



where d(j)k and 6(j)k sue Fourier transforms of the canon- 
ical field and of its conjugate momenta respectively, 
shifted by the "zero mode" 0o = {<P) and 0o — {4>)- This 
gives 



-I- id 



(B6) 



From the start we include in LJk the interaction with the 
bath of fluctuations. 



u;2 = fc2 + + + 3A(J02)^ 



(B7) 



i.e. Ofe correspond to "quasiparticles" . The second or- 
der correlators in homogeneous and isotropic background 
should be "diagonal" 



{alaq ) 

{akCLq) 



(B8) 
(B9) 



1. Microscopical Equations of Motion 

We derive equations of motion for the zero mode and 
for wave amplitudes starting from 

□0 -I- M^(t) + X(j)^ ^0. (BIO) 

a. Zero-Mode. Averaging IjBlOp we obtain 

4>o + + 3X{6(f>^))(Po + X4>1 + X{6(f>^) = . (Bll) 

In this equation (Scf)^) is small compared to the other 
terms and may be neglected locally in time in a state 

) is ei- 



(B3) which is close to a Gaussian. If additionally 



ther weakly varying or small compared to all other terms 
in (|Bllll the solution is given by the Jacobian cosine cn 



(po{t) ~ 00 cn fit, 



1 X<j>o 
V2 



where (j)o is the amplitude and 



fi = ^J X^l + 3X(5(j)^) + A'P . 



The period of this function is 



f 1 

V\/2 M 



(B12) 



(B13) 



(B14) 



where K{y) is the complete elliptic integral of the first 
kind. This defines the effective frequency, lUc — 27r/To. 
In the large amplitude limit, /i — A^^^^Oi we find uJc — 
For arbitrary one can write the following de- 
composition 



8 ^J, 



(B15) 



which is fairly accurate, the maximum deviation from 
the exact expression is 3% at /.t = A^/^c^q. For small 
amplitude of zero-mode oscillations, this expression can 
be further approximated as 



1 



3A0§ 



(B16) 



where M^s = {M^ + 3A(J02))i/2^ This deviates from 
exact result by less than 4% at X(f>o < M^g. 

For a general discussion of the kinetic equations in the 
background of a zero mode it might be useful to expand 
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4>o(t) in a Fourier series. However, this decomposition 
for the elhptic Jacobi function is strongly dominated by 
the first harmonic with uj = ujc- Even at /i = A^/^^O: the 
relative amplitude of the first harmonic is « 0.96, and it 
approaches unity with decreasing 0o- Therefore, in what 
follows we will restrict ourselves to the first term in the 
Fourier decomposition of (j)o{t). 

It is useful to define wave amplitudes for the zero-mode 



a, = V2^0oe-'"=* (B17) 
in terms of which the zero-mode can be represented as 

ar + a* 



(B18) 



One can also introduce an effective occupation number 
of " condensed waves" , 



2Wc<^n 



(B19) 



b. Wave amplitudes. The equations of motion for 
the wave amplitudes with non-zero momentum can be 
written as 



1 ^fe * , ^(3) 
Z ujk 



(4) 



(B20) 



where 



(3) _ 



„,)]<5W(fc-pi-P2) 



Ci^^ = -iX I dQki 
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IJpx^Vp2^ 



S^p,)]d^''\k-P,-P2-P3) 



d'^Pid'^P2 
/2^(27r)3'i/2 ' 



d'^pid'^P2d'^P3 
\/2ZJr(27r)5'i/2 



and 

(3) 

Cf, describes three wave interactions in the background 

of a zero-mode-and, while C^j^^ corresponds to four wave 
scattering. The averages in these expressions appeared 
because, first, we separated the zero-mode out of the 
equation for fluctuations, and, second, we used the ef- 
fective frequency for quasiparticles, Eq. HB7|) . Due to 
this choice the averages of C*^°^ times or will have 
the structure of cumulants, which in turn will deviate 
from zero only due to correlations induced by processes 
of scattering. 

Multiplying ljB20|l by or by au and subtracting the 
complex-conjugate expressions, we find 

flu = — Re cTfc -^Im/3(fc) -t-Im/4(fc) (B21) 

i&k = 2iOk(Jk + l:—nk+Il{k)+Il{k) (B22) 
2 uJk 



where 

h{k) = 6A0O j dnki2(al6(j)p,6<j)p.,)c (B23) 

h{k) ^ 2\j dnki23{al6cj)pMp2Hps)c (B24) 

Here (. . . )c denotes the cumulants, which in a diagram- 
matic language are identified with connected diagrams, 
see e.g. Ref. |93|. 

In situations when /3(fc) and Ii{k) are negligible, Eqs. 
(IB2ip and l|B22p describe particle creation in a time- 
dependent classical background, or parametric excitation 
if it is periodic. (Note that the quantum version of these 
equations at this stage can be obtained by the formal 
substitution — s- 1/2 -|- n^.) 

Note the following 

• In our case tuk contains an rapidly oscillating term, 
due to interaction with the zero mode. However, 
at late times and at large k it is small, e.g., in our 
numerical integration in the region of k near the 
peak of the spectral energy distribution, this term 
is of order 10"^, see Figs Q] and |21 We neglect this 
term in what follows. 



• The coefficient in front of the integral Eq. HB23|I 
is rapidly oscillating. Moreover, oscillations are 
not harmonic if the amplitude of 00 is large. Un- 
harmonicity can be dealt with by expanding (f>o{t) 
in Fourier time series and considering each of the 
terms separately. We restrict ourselves to the lead- 
ing harmonic in this expansion since at late times 
the unharmonicity is small. 

• The cumulants contain different combinations of a^ 
and a^, see Eqs. (IB4|) . It is well known that the 
leading contribution to the resulting kinetic equa- 
tion is due to the " resonant wave interactions" , 
or, in the language of particle physics, only those 
terms survive, which are on the "mass-shell" . In 
our case those will be {apap_^ap^)c and {apap^ap^)c 
for interactions which involve the zero mode, and 
{apap_^ap^ap^) c otherwise. We restrict our atten- 
tion to these cumulants only. 



• We neglect "anomalous" correlators, ak- These are 
small in the inertial range of turbulence as our lat- 
tice calculations show, but may be important oth- 
erwise. 



Leading Asymptotic of Collision Terms in 
Kinetic Approximation 



For a free random field the cumulants Eqs. I|B23|I and 
(|B24|) are zero, and hk = to the first order in perturba- 
tion theory. To calculate Uk in second order one has to 
know the solutions for cumulants in the first order with 
respect to interactions. 
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We use the equation of motion for wave amplitudes, 
Eq. HB20|I . to calculate the time derivatives of the cu- 
mulants, dta*{apap^apjc and i9t(a* ap^ap^OpJc. Higher 
order correlators which appear in this procedure can be 
used in zeroth order of perturbation theory, i.e. they can 
be decomposed in Uk assuming Gaussianity. To simplify 
the equations we use the following notations for products 



of Uk which appear in these decompositions 
^PiP2 = ricnp^np^ - UcUp {up^ + Up^) (B25) 

We find, keeping the terms which will have resonant be- 
havior 



dtaliapQp^apJc 
dt {a*pa*p^ap^ap^)c 



i{uJc + '^p-i^Pi-'^P2)a*c{a*pap^ap^)c + 



(B27) 



i6X (5('^)(p + pi - P2 - Pa) 



(B28) 



These equations have the common structure, iJ — IS.ljJ— 
A. Since A corresponds to the zeroth order in perturba- 
tions, it assumed to be time independent here. An Ap- 



propriate particular solution for cumulants is therefore 
given by J = A/(Atj -I- le), see e.g. Ref. [g^- Using the 
relation Im(a; -I- it)^^ — —'k5{x) we obtain 



(-27r)rf/2-i ^2cjc2tjp2tjpi 2u;p 



T */ * \ ^ o\ f2 r^y^c L^p Lup Lup /T3on^ 



Ini(apap.ap2«p3)c - 3A (2^)3./2-i^2c.,2c.,,2^,.2..,3 ^'^^ ^ ^ 



3. Isotropic Wave Kinetic Equations 



Applying this result to Eqs. HB23|I and IIB24|l we obtain 
the kinetic equation for wave occupation numbers Uk 



where 



7-(3) 

r(4) 



and 



dn'pp^ 



r(3) 



(B31) 



""piP2^ P1P2 



JQfc Pl pfcpi 
""P2P3 ^ P2P3 



F^p, (B32) 
(B33) 



_ 18A2 d<^pid'^p2 S'^'^Kk -Pl- P2) 
^ (27r)''-i2wc2wfc2wpi2cjp2 



18A2 d'^pid'^p2 S'-'^'^ (fc -f Pl - P2) 

{2T:Y~^2u;k2ujp,2ujp^2uJc 
X 6{ujk +ujpi - ujp2 - i^c) 



(B34) 



(B35) 



dfl'' 

P2P3 



_ 18A2 d'^pid'^p2d'% (5('^)(fe +P1-P2- Ps) 



{2T:Y<^-^2uJk2ujp,2uOp^2uOp^ 



X 5{ujk 



"P2 



- ^Ps) 



(B36) 



Both terms in HB31|I describe scattering processes of 
two waves into two other ones. In (|B32p one of them 
comes from the zero-mode, while in (|B33p all four are 
from the fluctuation field. Energy conservation in the 
interactions with the zero-mode, Wc+ti^fc— Wp^— — £k— 



be written as conservation of the energies 



(B37) 



which for small zero- mode amplitude (where ujc — i-^p=o) 
equals to the relativistic kinetic energy. Therefore, trans- 
port of energy over momentum space should be consid- 
ered as transport of kinetic energy in this case. 

These relations, Eqs. ||B32|| and ||B33)| . for 3- and 4- 
particle collision integrals can be reduced to one and two 
dimensional integrations respectively, if the distribution 
functions are isotropic, for details see e.g. Refs [gF 
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The collision integral for 3-particle interactions in d ■ 
3 takes the form 



r(3) 



d£2 - "fc ("3 + "2)] 



+ 2 



d£2 [Mnk + rii) -n,,n^] 



(B38) 



where Si = s{pi) and Ui = n{ei). Energy conservation in 
this case corresponds to £3 = Efc — £2 and £1 = £2 — £fc- 

The collision integral for 4-particle interactions in 
theory reduces to 



9A2 



32 TT"^ (^A; k 



poo poo 

/ d£2 / de^DF{n) , 
Jo Jo 



(B39) 



where D = mm[k,pi,p2,P3\ and £1 = £2 + £3 — £a; > in 
arguments of F{n) = (n^, + n^^) n2ft3 — nj^rii (77.2 + ^^.g). 

Note an interesting fact. Apart of the prefactor, 
Eqs. (|B39|) and (|B38|) are functions of relativistic kinetic 



energy, 7^; = {tokk) ^/(£fc). This gives for the flux of 
kinetic energy in 3 dimensions (cf. with Eq. H15|l^ 



d^p 



/(3) 



(27r)3 P P 
1 r^" 
'2^ 



(B40) 



deef{e) , 



where we have used pdp — code. Therefore, the turbulent 
flux should correspond to a particle distribution being 
a power law of relativistic kinetic energy. Remarkably, 
we have solved for the turbulent fluxes without the usual 
assumption of scale- independent dispersion law, Eq. (|17|l . 
In fact, the dispersion law was that of relativistic field 



theory, col = k"^ + M^f^. We do observe a single power 
law for the particle distributions as functions of kinetic 
energy in our lattice integration, even in situations when 
Mgff is large in the inertial range, see Fig. 01 the upper 
panel. 



4. Kinetic equation for zero-mode 

The kinetic equation for the wave occupation numbers 
has to be supplemented by the kinetic equation for the 
zero mode. We start with the equation 



60 + L^c'PO 



(B41) 



and repeat the procedure of the previous subsections. As 
analog of Eq. (|B21|) we obtain 



2 Aim 



(B42) 



Substituting (|B4I) and solving equation for higher order 
correlators we get 



d'^k 
(27r)° 



(3) 



(B43) 



This result is not surprising since 4-particle collisions 
conserve particle number. In the model we consider, 3- 
particle interactions are derived from the 4-particle col- 
lisions with one of the particles being replaced by the 
condensate. Therefore, Eq. Hij43|l can be interpreted as 
a conservation of the total occupation number, in parti- 
cles and in the condensate, Uc + (27r)~'^ / d'^k Uk = const. 
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